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ABSTRACT 

To better understand the impact of supernova (SN) explosions on the evolution of galaxies, 
we perform a suite of high-resolution (12 pc), zoom-in cosmological simulations of a Milky 
Way-like galaxy at 2 : = 3 with adaptive mesh refinement. We find that SN explosions can 
efficiently regulate star formation, leading to the stellar mass and metallicity consistent with 
the observed mass-metallicity relation and stellar mass-halo mass relation at 2 : ^ 3. This is 
achieved by making three important changes to the classical feedback scheme: i) the different 
phases of SN blast waves are modelled directly by injecting radial momentum expected at 
each stage, ii) the realistic time delay of SNe is required to disperse very dense gas before a 
runaway collapse sets in, and iii) a non-uniform density distribution of the interstellar medium 
(ISM) is taken into account below the computational grid scale for the cell in which a SN 
explodes. The simulated galaxy with the SN feedback model shows strong outflows, which 
carry approximately ten times larger mass than star formation rate, as well as smoothly rising 
circular velocity. Although the metallicity of the outfiow depends sensitively on the feedback 
model used, we find that the accretion rate and metallicity of the cold fiow around the virial 
radius is impervious to SN feedback. Our results suggest that understanding the structure of 
the turbulent ISM may be crucial to assess the role of SN and other feedback processes in 
galaxy formation theory. 
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1 INTRODUCTION 


Outflows are common in actively star-forming galaxies. Multi¬ 
wavelength images of a local starburst galaxy, M82, exhibit a 
clear bipolar outflow, with X-ray emitting gas at the centre (e.g. 
McCarthy , van Breugel & Heckman|ri987| [Strickland & Heck*^ 
m^2009| l. Lyman-break galaxies (LBGs) at high redshifts {z 


2 — 3) ubiquitously reveal significant blueshifted absorption of low- 
ionization metal lines, such as Cll A1334 ( jShapley et al.|2003^|Stei-| 
[del et al.|2010j ). These outflows often carry a large amount of gas 
that is comparable to or even several times larger than the mass of 
newly formed stars ( Martin|19^ [2005 1 jSteidel et al.|2010| [New-j 


[man et al.|2012[|Chisholm et al.|2014j ), indicating that stellar feed¬ 
back is the main process shaping galaxy properties. However, little 
is yet known about which types of stellar feedback mechanisms are 
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most important in regulating star formation in the baryon-driven 
picture of galaxy formation 

Studies show that a variety of galactic phenomena can natu¬ 
rally be explained by the presence of feedback from stars. jShetty &| 
|Ostrikerl ( |20121 ); [^m, Ostriker & Kiml ( l2013[ ) demonstrate that en¬ 
ergetic supernova (SN) explosions can drive turbulence in the inter¬ 
stellar medium (ISM), resulting in a low star formation efficiency of 
~ 1% per dynamical timescale of a galaxy, consistent with galactic 


scale observations (lKennicutt|1998||Bigiel et al. [2008 1 [Evans et al. 

2009 

1 . The evolution of the mass-metallicity relation {[Erb et al. 

2006 

IMaiolino et al.|20081|Zahid et al.|20131|Steidel et al.|2014J 


also seems to require that metals are not instantaneously recycled 
to form the next generation stars, but instead are first blown away 
from the galaxy ( |Dalcanton|2007[|Finlator & Dave|2008] l. Finally, 
semi-analytic calculations of galaxy formation suggest that, in or¬ 
der to match the galaxy luminosity functions below L* (e.g. [Bell 
et al.|2003 ► or the stellar mass to halo mass relation (e.g. [Most^ 
et al.f2010[ Guo et al.|2010j |, the regulation of star formation by SN 
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needs to be very effective so that only a small fraction ( ^ 20%) of 
baryons turns into stars (e.g. |Cole et al.|2000||Khochfar et al.|2007l 
ISomerville et al.|2008| references therein). 

Attempts have been made to understand the formation of 
galaxies with SN feedback in a fully cosmological context using 
numerical methods ( |Katz||199'^ |Cen & Ostriker|p^^ [Navarro] 
|& White||1993| ). However, [Navarro & Steinm^ ( |2000| ) find that 
simulated galaxies with kinetic feedback rotate too fast, inconsis¬ 
tent with the local Tully-Fisher relation { [Tully & Fisher[[TW7] l. 
The angular momentum catastrophe is attributed in part to the fact 
that dynamical friction of infalling gas is overestimated as stan¬ 
dard smooth particle hydrodynamics (SPH) cannot capture hydro- 
dynamic instabilities at contact discontinuities accurately and gas 
comes in cold and clumpy ( [Navarro & Steinmetz[[2000l [Sijacki] 
[et al.|2012] ). However, simulations with grid-based codes and weak 
SN feedback also find that galaxi es are too compact ([Joung, Cen &| 
[Bryan|200^ [Kimm et al.|20fT^ [Hummels & Bryan|2012| l, as gas 

cools excessively and forms stars at the galaxy centre. How to over¬ 
come the overcooling problem ( |Katz|199^ and robustly include the 
effects of SNe are important tasks for all numerical galaxy forma¬ 
tion studies. 

In order to ensure that energy from SNe is not radiated 
away artificially, several authors introduced the cooling suppres¬ 
sion model in which gas near young stars is assume to be adia¬ 
batic for several to tens of Myr ( jGovernato et al.[[2010| [Guedes] 
[et al.|20lT] l. With this feedback model, [Guedes et al.| ( |2011| l find 
that a realistic disc galaxy with rising rotation curves can be pro¬ 
duced within a ACDM paradigm if stars are permitted to form 
only in dense environments (nn = 5 cm“^), so that star formation 
and outflows are strongly clustered. They claim that the collective 
winds are efficient at removing low-angular momentum gas that 
is accreted at high redshifts (see also [Brook et al.|2()TT] l. [Oh et al. 

( 2011 1 analyse the simulation of a dwarf galaxy by |Governato et^ 
( 2010 i, and find that the simulated inner slope of mass density as 
well as the shape of rotation curves are consistent with those of 
THINGS galaxies ( [Walter et al.pOOS] ), demonstrating that strong 
feedback is crucial to reproduce the observations. Yet the timescale 
imposed for the adiabatic phase is controversial, given that it relies 
on the density of the ISM, which is largely under-resolved in the 
simulations. 

Recently, other forms of stellar feedback have been suggested 
as important mechanisms to control the growth of galaxies. [Walch] 
[et al.[ ( [20T^ ; [Dale et al.[ ( |2014j l show that photoionisation is capable 
of dispersing a small cloud of mass 10^ M© in ~ 3 Myr by gener¬ 
ating an over-pressurised HII region (see also [ Vazquez-Semadeni| 
[et al.|20I(}}[Sales et al.|2014j l. [Dale et al.] ( [2014j ) point out that winds 
from massive stars can create a cavity of ~ 10 pc at the centre of 
10^-10^ M© clouds, but they alone cannot unbind the gas from the 
system (see [Geen et al.[2015] for a uniform case). One also expects 
that radiation pressure from ionizing photons has a similar effect 
as stellar winds (for metal-rich gas), given that their momentum in¬ 
jection rates are comparable ( [Leitherer et al.[199^ . However, [Sales] 
[et al.[ ( [20T4] ) assert that the inclusion of radiation pressure by ion¬ 
ising photons on top of photoionisation has a subdominant effect, 
because gas near a massive star is transparent to the Lyman-limit 
photons once ionised and cannot be accelerated continuously (c.f. 
[Wise et al.[2012b[ ). The presence of dust can change the conclusion 
though, as it can also absorb the ionising radiation. More impor¬ 
tantly, when gas is optically thick to infrared (IR) photons, the lat¬ 
ter can in principle be multiply scattered supplying 10-100 times 
the original momentum of the ionizing radiation ( [Murray, Quataert] 
[& Thompson[2010[ ). Hopkins and collaborators conduct systematic 


studies to gauge the relative importance of photo-heating, radiation 
pressure, and SN explosions in isolated discs ( [Hopkins, Quataert"^ 
[Murray[2011 1 [2012a|b[ l, and find that radiation pressure with a mo¬ 
mentum boost of 10-100 is required to drive violent galactic out- 
fiows in actively star-forming galaxies (c.f. Rosdahl et al.|2015|l. A 
similar conclusion is reached by|Agertz et al.|(|2013|l;|Aumer et al.[ 

( |20T^ . 

It should be noted, however, that several inaccurate assump¬ 
tions are made concerning SN explosions in previous cosmologi¬ 
cal simulations. First, even though the initial expl osion energy of 


Nomoto et al. 


10 erg is used, as guided by observations (e.g. _ 

[1993| ), the momentum transfer from SN to the ISM is likely 
to be underestimated due to limited resolution. Because a large 
amount of energy is deposited in a very concentrated form, SN 
blast waves must undergo an adiabatic phase during which ra¬ 
dial momentum should increase by up to an order of magnitude 
(?^ 3 X 10^kms“^M©), compared to the initial ejecta momen- 


tum (|Sedov|1959| Taylor|1950l Chevalier|1974[[Cioffi, McKee &[ 


[B ertschinger [[ 1988 [Blondin et al.[[1998| l. Without modelling this 

phase properly, hydrodynamic simulations underestimate the mo¬ 
mentum transfer from SN. [Thornton et al.] ( [1998[ ) performed a pa¬ 
rameter study of the explosion in uniform media using one di¬ 
mensional hydrodynamic simulations, and showed that the final 
radial momentum is a weak function of the bac kground density 
(oc and metallicity (oc (see also|Kim & Ostriker 


[2015[|Geen et al.|2015[ ). The final momentum transferred to an in¬ 


homogeneous medium turns out to be surprisingly similar (Kim 


& Ostriker [2015[ [iffrig & Hennebelle[[2015[ [Martizzi, Faucher- 


Giguere & Quataert[2015[ l, as it has a weak dependence on density. 
Although this has been appreciated in studies of the generation of 
turbulence in the ISM (e.g. [Shetty & Ostriker|2008[ l, only recently 
a few cosmological simulations begun to take into account the mo¬ 
mentum at different phases of SN expansions ( [Kimm & Cen|2014| 
[Hopkins et al.[2014] ). Second, each star particle is often assumed to 
release SN energy in a single event at a fixed time delay of ~ 10 
Myr ( [Dubois & Teyssier|2008[ l. This may overestimate the impact 
of SN, as clustered explosions can work together to create galac¬ 
tic fountains. But at the same time, it is likely that the regulation 
of star formation is mistimed when the delay for SN explosions 
is fixed to a single time such as 10 Myr. [Aumer et al.[ ( [2013[ ) re¬ 
ported that in the absence of feedback mechanisms active before 
SNe, to regulate star formation in their simulations, they needed to 
use an early time delay of 3 Myr. This time delay appeared to leave 
galaxy star formation histories unchanged when radiation pressure 
was included in the simulations. Lastly, even though up-to-date 
cosmological simulations can now employ a spatial resolution as 
high as 10 parsecs, this is not enough to resolve the turbulent struc¬ 
ture of the ISM where the volume-filling density is smaller than the 
mean (e.g. [Federrath & Klessen|201^ . [Iffrig & Hennebelle[ ( |2015] ) 
demonstrate that the expansion of a SN out to some fixed radius 
depends on the local environments of the explosion. This implies 
that determining the mass swept up from the host cell of a SN may 
not be trivial. In this respect, questions remains regarding how mo¬ 
mentum should be injected on resolved scales. 

The aim of this paper is to understand the impact of SN by 
realistically modelling the momentum injection from SN. Can the 
momentum input from SN alone possibly explain the observed re¬ 
lations, such as the mass-metallicity relation or stellar mass to halo 
mass relation? Why do the classical feedback schemes fail to repro¬ 
duce these? In order to address these questions, we perform a suite 
of zoom-in, cosmological hydrodynamic simulations with adaptive 
mesh refinement (AMR). We describe our simulation set up and a 
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High-z galaxies with mechanical supernova feedback 3 


variety of feedback models tested in this work in Section 2. Our 
results on the evolution of stellar mass, metallicity, and kinematic 
properties are presented in Section 3. We discuss the effects of each 
prescription and why some feedback models are more efficient at 
suppressing star formation in Section 4. Finally, we summarise our 
findings in Section 5. 


2 SIMULATIONS 


We use the tree-based Eulerian AMR hydrodynamics code, 
Ramses < |Teyssier||2002| ), to study the impact of SN feedback on 
the evolution of a Milky Way-type galaxy. The initial conditions 
for the simulations are identical to those used for the Nut suite 
presente d in [Powell, Slyz & Devriendt| ( |2011| ); |Geen, Slyz & De-| 
|vriendt| (|2013t, which were generated using the MPGRAFIC soft- 
ware ( [Prunet et al.|2008||Bertschinger|2001| l with WMAP5 cosmo¬ 


logical parameters (Qrn = 0.258, Qa = 0.742, Qb = 0.045, erg = 
0.8, iTo = 72kms“^ Mpc“^) jPunkley et al.||200^ . The com¬ 
putational domain represents a periodic box of 9 h~^ Mpe on a 
side, covered with 128^ cells. Within this volume, three nested 
grids are placed in a spherical region of radius 2 Mpe (comoving) 
encompassing a dark matter halo of Mvir = 5.5 x 10^^ M© at 
z = 0. The corresponding mass resolution of a dark matter particle 
is 5.5 X 10^ M©. Only the spherical region is allowed to be further 
refined. Figure shows the composite images of density, tempera¬ 
ture, and metallicity around the Nut halo at three different redshifts 
{z =3, 5, and 7) from a run with SN feedback. 

The Poisson equation is solved using a multi-grid method 
( |Guillet & Teyssier||2011| ) for grids with level 7 < Iv <9, 
while a conjugate gradient method is used for more refined grids 
{Iv > 10, |Teyssier|[2002| ). The Euler equations are integrated in 
super-comoving coordinates ( Martel & Shapiro] 1998| ) with a HLLC 
Riemann solver ( |Toro, Spruce & Speares|1994| ). We adopt the typ¬ 
ical courant factor of 0.7. 

Metal- dependent radiative cooling is modeled with cooling 


curves from Sutherland & Dopita |l993| ), down to ~ 10^ K. Gas 
cooler than 10^ K can lose energy further via metal fine struc¬ 
ture transitions, following |Rosen & Bregman] { |1995| ). A uniform 
ultraviolet background is turned on at z = 8.5 using [Haardt &| 
|Madau| ( p996| ), and modeled as a heating term in the energy equa¬ 
tion. However, gas denser than nn = 0.01 cm“^ is assumed to be 
self-shielded from the radiation ( [Rosdahl & Blaizot|2012|[Faucher-| 
[Giguere et al.|20l6] ). The simulations start with gas of metallicity 
Z — 2 X 10“^ at z = 499, as the primordial SN explosions in 
mini-haloes are not resolved in our simulations (e.g. | Whalen et al.| 
iTOOSlIMsTet al.|2012b| l. 


In order to better resolve the ISM in galaxies, we employ a 
maximum physical spatial resolution of 12 pcQ(= 12.5 Mpc/2^°) 
by triggering 10 more levels of refinement when baryonic plus dark 
matter mass in each cell becomes greater than 4.4 x 10^ M©. In 
addition, we also enforce that gas denser than nn = 1, 42.5, 
340 cm“^ is always resolved at least to resolutions of 48 pc, 24 
pc, and 12 pc cell, respectively. At the final redshift {z = 3), the 
main halo is resolved with 3.0-3.7 x 10® computational cells, de¬ 
pending on the feedback model. We find that 59%, 14%, and 5% of 


^ We perform a convergence test by comparing the MFBmp run and its 
twin version with a higher resolution (6 pc) down to 2 ; = 5.5, while keeping 
other parameters fixed, except for the mass of a star particle (203 M©). We 
find that the stellar mass of the main galaxy is converged within 10% at 
z = 5.5. 


the gas mass inside 0.1 i^vir (~ 3.6 kpc ~ 2Re) is resolved on 48, 
24, and 12 pc resolution elements in the case of the NutMFBmp 
run. The main galaxy in the NutTFB run uses more high-resolution 
cells (40%, 32%, and 24% of the gas mass is resolved on at least 
48, 24, and 12 pc cells, respectively), as the ISM becomes denser 
than that of MFBmp due to the overcooling. 

Dark matter haloes are identified with the HaloMaker { |Tweed| 
|et al.|200>^ . The virial radius of each halo is defined as the radius 
within which the mean density is equal to the critical density of 
the Universe (perit = 377(2:)^/8 7rG) times the virial over -density 
(Acrit = 187r^ + 82x — 39x^, Bryan & Norman|ri998| l, where 
X = dm/(dm + u^dA) — 1 and a is the scale factor. 


2.1 Star formation 

We model star formation as a stochastic process (e.g. |Katz 
[Rasera & Teyssier|2006| l, based on a Schmidt law < |Schmidt|1959| ), 

Pstar — CffPgas/^ffj (1) 

where pgas is the gas density, tff = ^37r/32Gpgas is the free-fall 
time of the gas, and Cff is the efficiency of star formation per free- 
fall time. Specifically, we examine the probability of forming a star 
particle in a cell denser than a certain threshold density (ttsf) in 
each fine time step. We adopt nsF = 42.5 cm“®, motivated by the 
Larson-Penston density (plp — 8.86 c ^/tt G Sx‘^) above which 
isothermal gas undegoes a runaway collapse ( |Larson|1969||Penston| 
|1969[[Gong & Ostriker|2011||2013| ), where Cs is the sound speed 
and (5a; is the scale of the gas cloud. For a typical temperature of the 
cold ISM (~ 30 K) and the finest size of our computational grid (12 
pc), this yields nn = 42.5 cm“®. Although the threshold density 
is slightly lower than the typical mean density of a giant molecular 
cloud (nn ~ 100 cm“®), we note that the majority ( ^ 70%) of 
stars form in such environments, as star formation is generally a 
slow process ( |Kennicutt|1998[|^umholz & Tan|2007[[Evans et al.| 
|2009| see below). 

The number of stars formed (Np) is determined by drawing a 
random number from a Poisson distribution with a mean of 

^ = eff(pAx /m^^rnin) (Atsim/tff) , (2) 

where Ax is the size of a computational cell, and ?n*,min = 
a Np nsF mu Ax^^^/Xu is the minimum mass of a star particle. 
Here mu is the hydrogen mass, Ah (=0.76) is the mass fraction of 
hydrogen for the primordial composition, and a is a. parameter that 
controls the minimum mass, which we use a = 0.264 to achieve 
^★,min = 610M©. We assume that 2 per cent (eff = 0.02) of 
the dense gas (nu >42.5 cm“^) turns into stars per free-fall time 
UKrumholz & Tan|2007| l. 


2.2 Feedback models 

We examine six different feedback models to investigate the im¬ 
pact of core-collapse SN on the evolution of galaxies. We use the 
Chabrier initial mass function (IMF, |Chabrier et al.|2005| l through¬ 
out this study. We assume that 31 % of the total mass is lost via Type 
II SN explosion. For a progenitor mass of 19.1 M©, this means that 
a star particle of 610 M© would host 10 SNe. The corresponding 
specific energy per M© at the time of explosion is 1.6 x 10^® erg. 
The metallicity of the ejecta is taken to be 0.05, which is equivalent 
to a yield of 0.016 for a simple stellar population with 1 M©. 

Feedback models may be distinguished by how the energy or 
momentum is put into the ISM. Note that we do not turn off the 
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4 T. Kimm et al. 



Figure 1. Composite images of the Nut galaxy at 2 ; = 7, 5, and 3 from the NutMFBmp run. The pink, blue, and green colours represent the projected 
distributions of density, temperature, and metallicity, respectively. The image measures 2.715 Mpc (comoving) on a side. The white circle denotes the virial 
radius of the host halo of the Nut galaxy (41 kpc at 2 : = 3, physical), which is fed by roughly three large-scale filaments. It can be seen that the metal-enriched 
SN-driven winds extend out to several virial radii of the halo at 2 ; = 3. 


radiative cooling ( |Mori et al.|1997[[Stinson et al.|2006j) or hydrody- 
namically decouple winds from the ISM ( [Scannapieco et al.|2006[ 
|Oppenheimer & Dave|2006| l. We describe the details of each model 
below. 

(i) Reference run (NutCO) 

The reference run (NutCO) does not include any energy or momen¬ 
tum feedback, but gas and newly processed metals from SNe are 
added to their host cell, once a star particle becomes older than 10 
Myr (so called “metal feedback”). Since there is no direct energy 
input in this case, star formation is essentially controlled by the cos¬ 
mic inflow and uniform background UV heating. Note that we use 
a higher minimum mass for star particles (?Ti*,min = 2310 M©) 
than other runs to reduce the computational cost. 

(i) Thermal feedback (NutTFB) 

The total energy of 10^^ erg per SN is simply added to the host cell 
of the star particle, along with ejected gas and metals. Each star 
particle is assume to generate a single SN event at 10 Myr after the 
birth. This can increase the temperature of a cell up to ~ 3 X 10® K 
in a low-density medium (i.e. Mgas ^ Mgj). In principle, the over¬ 
pressurised gas can drive the expanding motion depositing ~ 30% 
of the total into the kinetic energy ( |Chevalier|1974| ). However, as is 
well known in the literature, the atomic and metal cooling processes 
can rapidly radiate the internal energy away before the SN-driven 
blast wave sweeps up the surrounding medium in dense environ- 


ments or in simulations where the cooling radius is under-resolved 

^Katz| 

119921 [Navarro & White|1993| Abadi et al.|2003|[Slyz et al.[ 

2005 

[Hummels & Bryan|2012|[Kimm & Cen|2014[>. To circum- 


vent this issue, some studies assume that gas is adiabatic for several 
tens of Myr (|Mori et al.| 19^[Thacker & Couchman|200T||Stinson| 
|et al. 12006) |Governato et al.|2007[|Teyssier et al.|2013| l. Our ther¬ 

mal feedback model does not include the prescription to properly 
address excessive radiative losses. 

(ii) Kinetic feedback with a small mass-loading (NutKFB) 

We also test the kinetic feedback scheme of [Dubois & Teyssier| 
^2008t , in which the initial explosion energy of E^sn — 10^^ erg 
is modelled by increasing the kinetic energ30 of the surrounding 

^ When the momentum added is cancelled out, the kinetic energy is used 


gas. An important parameter in this model is a mass-loading factor 
{rjw = Mioad/Mej), which accounts for how much gas is entrained 
with the SN ejecta. Notice that a mass-loading factor (77) is often 
defined as the mass ratio between the outfiow rate (Mout) and the 
star formation rate (Mstar), which is different from what we use 
in this study. Empirical determinations of the mass-loading factor 


suggest that it varies depending on the mass of a galaxy 

Martin 

20051 

[Rupke, Veilleux & Sanders[2005|[Weiner et al.[2009 

Martin 

et al.| 

I2OI2I [Newman et al.[[2012l [Chisholm et al.[[20141 

Arribas 

et al.[[2014||, but is roughly on the order of pw ~ 10. Motivated 


by this, we adopt r]w,ma^ = 10 (i.e. 77 ~ 3). Notice, however, 
that r]w can be as small as zero if there is little gas in the host cell 
of SN. The velocity of the neighboring gas is then determined by 
V = ^/2 E^sn/(1 + Vw) A7ej, hence the input momentum can be 
as 3.3 times as high as the initial momentum if there is enough gas 
to entrain from the host cell. Again, SNe are modeled as a single 
event at 10 Myr per individual star particle, but this is evaluated at 
every coarse time step (~ 0.16 Myr) to reduce the message passing 
interface communications between computing cores. 

(iii) Mechanical feedback (NutMEB) 

[Kimm & Cenj ( |2014j ) introduced a physically based feedback 
model, which can properly explain the transfer of momentum at all 
stages of a Sedov-Taylor blast wave (from the adiabatic to snow¬ 
plough phase) (see also [Hopkins et al.|2014| l. The main difference 
from the kinetic feedback model is that the input momentum is cal¬ 
culated according to the stage of the Sedov-Taylor blast wave. Al¬ 
though the momentum available during the free expansion phase is 
only 4.5 x 10^ km M©, it increases as more gas is swept up 
by the shock (p oc ^/m). When the energy lost at the shell becomes 
significant (AE/E ^ 0.3) and the pressure inside the bubble be¬ 
comes comparable to that of the pre-shock region, the expansion of 
the shell becomes momentum-conserving. At this point, the shell 
mass gain is much higher than the decrease in velocity resulting in 
the net increase of the total momentum by a factor of ~7-10 at the 
end of the adiabatic phase. This factor can be smaller if the number 
of SNe is larger or the background density of gas is higher or the 


to increase the temperature to conserve the total energy. This is also the case 
for the mechanical feedback runs. 
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metallicity is higher, as the radiative cooling becomes more signif¬ 
icant. 

Specifically, the input momentum during the snowplough phase 
is calculated as ( [Blondin et al.|1998] Thornton et al.|199^ see also 
|Kim & Ostriker|2()T5 JGeen et al.|2015]) 

PSN,snow « 3 X 10®kms-^ Mo (3) 


where E^i is the number of SNe, nn is the hydrogen number den¬ 
sity, and Z' is the metallicity in solar units, with a minimum of 0.01 
(Z' = max [ZjZQ^ 0.01]). Whether the expansion is in the adia¬ 
batic or momentum-conserving phase is determined by the mass 
swept up from the host and neighboring cell. This is done by cal¬ 
culating the mass ratio (x) for each 48 neighbours (= A^nbor) per 
explosion as 


X — C^Afswept / , 


(4) 


where 


dMej = (1 - ^sn)Mej/iVnbor, (5) 

jfl/T _ f 1 (1 -/?sn)pho3tAa:® 

dAfswept — Pnbor I 2 ) ]\[ , (6) 

Ax and phost is the size and gas density of the host cell, respec¬ 
tively, and Anbor = 48 is the number of neighbouring cells (see 
Figure 15 in |Kimm & Cen||2014| l. Here /3sn determines the mass 
fraction of (Mej + phostAx^) to be left in the host cell. We take 
/3sn = 4/52 to distribute the mass as evenly as possible to the host 
and neighboring cells when the cells are on the same level. If x is 
greater than the transition mass ratio, 

Xtr = 69.58 (7) 

we input psN,snow (Equation [^, whereas momentum appropriate 
to the adiabatic phase is injected otherwise, 

_ J PSN,ad = V^2x Afej fe F/SN (x ^ Xtr) 

PSN \ / \ \ 5 v®/ 

( PSN,snow (X _ Xtr) 

where /e = 1 — is used to smoothly connect the two 

regimes. Note that the local properties of gas at which a SN ex¬ 
plodes can be significantly different from those of gas at which the 
star particle is formed, as we allow for the ISM to dynamically 
evolve for 10 Myr before the SN explosion. 


(iv) Mechanical feedback with multiple explosions (NutMFBm) 
This model is based on the same mechanical feedback scheme as 
NutMFB, but differs in that SNe from a star particle are distributed 
in time, as in the real Universe. In cosmological simulations with 
finite resolution, an individual star particle represents a coeval star 
cluster, which is usually assumed to have a single explosion at some 
fixed age. Our models (NutCO, TFB, KFB, MFB) also assume that 
SNe explode at 10 Myr, meaning that tens of SNe explode simulta¬ 
neously. In reality, some SN explosions commence at as early as 3 
Myr, whereas an 8 M© star evolves off the main sequence only at 
~ 40 Myr. The multiple SNe model is designed to mimic the dis¬ 
crete nature of the multiple explosions by randomly sampling the 
lifetime of massive stars for each particle. To do so, we use the in¬ 
verse sampling method for a polynomial fit to the integrated rate of 
SN Type II (Figureleft panel), 

n(< t) = vj T ^6 + c(.Z) (9) 

where te = t/lO^r, a{Z) = -0.268 + 0.139 Z' - 0.575 
h{Z) = 0.042 + 0.022 Z' + 0.0790 Z'^ c(Z) = -0.086 - 
0.170 Z' + 0.187 Z'^ and Z' = max(Z, 0.008)/0.02. We use 




t [Myr] t [Myr] 

Figure 2. Left: Integrated SN rates based on the population synthesis code 
Starburst99 (solid) and the corresponding polynomial fits (dashed lines, 
Equation [^. Different colours indicate different metallicities, as shown in 
the legend. Right: An example of the random sampling of discrete SN ex¬ 
plosions for a 1000 Mq star particle. 


the population synthesis code Starburst99 ( [Leitherer et al.|199'^ to 
compute the SN rates as a function of time and metallicity. Figure]^ 
(right panel) shows an example of the random sampling and corre¬ 
sponding mass evolution of a 1000 M© star particle. Note that, 
unlike NutKFB and MFB, we evaluate the explosion at every fine 
time step (several hundreds to thousand years) for a better accuracy. 


(v) Mechanical feedback with multiple explosions in a porous 

ISM (NutMFBmp) 

Particular attention is paid to how momentum is distributed to the 
neighboring cells in this model. In NutMFB and NutMFBm, all of 
the gas in the host cell of SN is re-distributed uniformly to its sur¬ 
roundings. This may be a reasonable assumption for an ISM with a 
smooth structure or for an ISM with a low Mach number where the 
volume filling density is not very different from the average den¬ 
sity ( |Vazquez-Semadeni|1994[|Federrath & Klessen|2012|l. How- 
ever, high-z galaxies are normally turbulent (e.g. [Forster Schreiber| 
let al.|2009| ), and our resolution elements may not fully resolve the 
porous structure of the turbulent ISM. If a SN explodes in a porous 
medium, energy will preferentially propagate through low-density 
channels ( Iffrig & Hennebelle|2015||Kim & Ostriker |2015||Mar-| 


[tizzi, Faucher-Giguere & Quataert|2015| t, unlike SN bubble expan- 

sion in a uniform medium. In this regard, it is likely that too much 
gas is entrained from the host cell of the SN in simulations with lim¬ 
ited resolution (including ours), leading to the underestimation of 
the initial gas outflow velocity in the neighboring cells. For exam¬ 
ple, if a SN explodes in a cell with nn = 10 cm“^ resolved to 24 
pc and a neighbouring cell has a density of 1 cm“^, the initial wind 
velocity would be about 30kms“^ for a uniform ISM. However, 
if the ISM were highly turbulent with a Mach number of 10, as 
is often observed in the Galactic star-forming clouds (e.g. |Roman-| 
|Duval et al.|2010| |, the volume filling density would be roughly ten 
times smaller for a collapsing cloud { [Molina et al.|2012[|Federrath| 
|& Klessen|2012| > and the wind velocity would therefore be roughly 
two times faster than the estimate with the uniform case. Note that 
the increase in the wind velocity is partially attributed to the fact 
that the input momentum would be 30% larger in the porous case. 
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Table 1. Summary of nut42 simulation parameters and physical ingredients. From left to right, columns are as follows: simulation name, type of SN 
feedback, minimum grid size (physical units), mass of a dark matter particle, minimum mass of a star particle, threshold density for star formation, star 
formation efficiency per free-fall time, final redshift of each simulation, time delay of a core-collapse SN, and remarks. ? 7 t(;,max is the maximum mass-loading 
factor with respect to the ejecta mass, used in |Dubois & Teyssi^|2008| DT08). /w,host indicates the fraction of the mass entrained from the host cell of a SN. 


Simulations 

SN II feedback 

Aa^rxiin 

rUDM 

'mstar,min 

nth 

eff 

^end 

fdelay 

Remarks 



[pc, physical] 

[M©] 

[M©] 

[H/cm3] 



[Myr] 


NutCO 

metals only 

12 

5.5 X 10^ 

2310 

42 

0.02 

3 

10 


NutTFB 

thermal 

12 

5.5 X 10"^ 

610 

42 

0.02 

3 

10 


NutKFB 

kinetic (DT08) 

12 

5.5 X 10^ 

610 

42 

0.02 

3 

10 

ryt(;,max — 10 

NutMFB 

mechanical 

12 

5.5 X 10"^ 

610 

42 

0.02 

3 

10 


NutMFBm 

mechanical 

12 

5.5 X 10^ 

610 

42 

0.02 

3 

realistic 


NutMFBmp 

mechanical 

12 

5.5 X 10"^ 

610 

42 

0.02 

3 

realistic 

porous ISM (/w,host = 0.1) 


as the background density is lower. In other words, under-resolving 
the host cell of SN can artificially smooth out the winds, reducing 
the initial velocity of the neighboring gas, especially along low- 
density channels. 

In an attempt to examine the effect of the porous ISM and esti¬ 
mate its possible consequence on the impact of SNe, we allow for a 
faster wind by reducing the amount of gas mass entrained from the 
SN cell. We simply assume that only 10% of the gas from the host 
cell of SN (/w,host = 0.1) is directly involved in the expansion. 
This means that Equation [^reads 

dMswept = Pnbor (^) ^ + + dMej. ( 10 ) 

We emphasise that this number is rather arbitrarily chosen, hence 
the comparison with observations should be taken with caution. 
Nevertheless, we note that it can clearly demonstrate the impor¬ 
tance of detailed modelling of momentum injection and address 
the question as to whether the momentum budget from SN suffices 
to regulate star formation. 

The parameters used in the six runs are summarised in Table[^ 


3 RESULTS 

In this section, we compare the baryonic properties of simulated 
galaxies to evaluate the effect of different SN feedback models. 
Our zoom-in region contains 99 dark matter haloe^with Mvir > 
10^ M© at z=3, including the progenitor of a Milky Way-like 
galaxy. Our analysis mostly focuses on the most massive galaxy, 
but we also present results for other smaller galaxies to see a gen¬ 
eral trend of the physical properties as a function of galaxy mass. 

The centre of a galaxy is often taken as the densest region 
of a dark matter halo in large-scale simulations (e.g. [Tweed et al.| 
|2009| ). But it is not necessarily the case in high-resolution simula¬ 
tions where mergers are more frequent and a galaxy is often very 
clumpy. Starting from the stellar plus gas distributions within 0.5 
i^vir, we redefine the galaxy centre by iteratively calculating the 
centre of baryonic mass within its half-mass radius (reff,m) un¬ 
til converged. This ensures that the centre of a galaxy is normally 
taken as the centre of the most massive clump. The galaxy stellar 
mass is then computed by summing up the mass of the star particles 
within 0.2 i^vir from the centre. 

The main galaxy in our simulation is hosted by a dark m at- 
ter halo with 5.5 x 10^^ M© at z = 0 ( [Kimm et al.||201Ia| ). At 


^ Subhaloes are not included in this sample. 



Figure 3. Halo merger tree of the Nut galaxy. The leftmost branch corre¬ 
sponds to the merger history of the main galaxy. The digits shown in the 
right side of a merger event indicate a merger ratio (7 = Mcen/^sat) be¬ 
low 20. The ratio for the very minor merger (20 < 7 < 100) is omitted. 
Major (7 < 4) and minor mergers are displayed as thick and thin solid 
lines, respectively. A bigger symbol indicates a more massive halo. Note 
that actual galaxy mergers take place later than the halo merger. 


2: = 3, its virial mass is 10^^ M©, which is the typical host 
halo mass of Lyman alpha emitters (LAE, [Gawiser et al.|2007j l. As 
a comparison, LBGs are thought to reside in more massive haloes 
with Mhaio ^ |Giavalisco & Dickinson|2001 Arnouts et al. 


|2002[[Melberger et al. [2005 [[Lee et al.|2009 ^elby et al.|2013j) 

Using the two Millennium simulations (jSpringel et al. 12005] 

[Boylan-Kolchin et al.||2009| l, [Fakhouri, Ma & Boylan-Kolchinj 
( |2010| | measured that a halo of mass 10^^ — 10^^ M© would typ- 
ically have ~2, 10, and 20 mergers of the ratio 1:3.3, 1:33, and, 
1:100 at 3 < z < 9, respectively. The Nut halo undergoes 6 mi¬ 
nor mergers (mass ratio of 4-20), and 15 very minor (mass ratio 
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MFBmV' 


Figure 4. Projected distributions of the density (first and third columns), metallicity (second column), stellar density (fourth column), and star formation rate 
density (fifth column) of the Nut galaxy in runs with different feedback models, as indicated in the top left corner of each panel. The dotted circles in the first 
two columns denote the virial radius of the dark matter halo. The last three columns show the central region of the halo. Note that all images are produced along 
the same direction (y-axis) for direct comparisons. The stellar component is markedly different from each other depending on the feedback model adopted. 
The spin direction of the galaxy is particularly sensitive to the choice of the feedback model. The strong outflow in the NutMFBmp run prevents the formation 
of a well-ordered gaseous disc at 2 : = 3. 
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z=3 4 5 6 10 15 



z=3 4 5 6 10 15 



Figure 5. Comparison of the galaxy stellar mass (left) and star formation history (right) of the Nut galaxy. The abscissa indicates the age of the universe 
(til). Different colour-codings correspond to different feedback models, as indicated in the legend. The right panel shows star formation histories calculated 
using star particles located within 0.2 virial radius of a dark matter halo at 2 : = 3. Note that this includes stars formed in all progenitors. It can be seen that 
mechanical SN feedback is more effective at suppressing star formation than either thermal or kinetic feedback. Star formation is very episodic in the MFBmp 
run, as SNe efficiently destroy star-forming clouds. 


of 20-100) mergers at 3 < z < 10, hence it has a relatively quiet 
merger hist ory (Figure [^. Not e that the only major merger occurs 
at z = 2.2 ( |Kimm et al!j201 la| ). In Figure]^ we show the projected 
distributions of gas density, metallicity, stellar density, and star for¬ 
mation rate density of the Nut halo and galaxy at z = 3 in each 
run. 


3.1 Star formation histories 


We begin our investigation of the evolution of galaxy properties 
with the cooling run (NutCO) in which energy from SNe is ne¬ 
glected. Because there is no pressure support from stars and the UV 
background radiation is self-shielded above 0.01 H/cc ( |Faucher-| 
[Giguere et al.||2010} |Rosdahl & Blaizot||2012| l, star-forming gas 
clumps experience a runaway collapse converting a large fraction of 
the total baryons into stars at all redshifts (Figurethe left panel). 
The growth of the stellar component closely follows the assembly 
history of its host dark matter halo (shown as a dashed line), in¬ 
dicating that star formation is limited by gas accretion. The main 
progenitor of the Nut galaxy forms stars at a rate of ~ 10 M© yr“^ 
in the last ~ 1.5 Gyr (Figurej^ right panel), resulting in a final stel¬ 
lar mass of Mstar = 1.9 X 10^° M© at z = 3. As a consequence 
of the efficient star formation, the galaxy ends up being gas-poor 
at z = 3 (/gas — Afgas/ [Afgas Afstar] ~ 0.15), inconsistcnt 
with the observations at z ~ 2 — 3 that gal axies in a similar mass 
range are typically gas rich (/gas ^0.5) (e.g. 

ITacconi et al.|20 !()] ). 

re|^s 


Mannucci et al. 2009 


Figure^shows the mass of the progenitor dark matter halos in 
which each star particle is formed. The plot demonstrates that most 
of the stars populating the Nut galaxy z = 3 in the NutCO run 
are formed in massive haloes. The amount of stars formed between 
3 < z < 6 (when the halo mass is 10.5 < logMhaio/M© < 
11 ) accounts for half of the total stellar mass, which refiects that 
more massive galaxies show a higher level of star formation activity 


( |Salim et W. |2007[ |Elbaz et ^|2007[ [Noeske et al.||2007 |Daddi 


et al.||2007l [Gonzalez et al.||20T0l [Whitaker et al.|[20T^ [Salmon 


et al.|2015] l. Thus, feedback must regulate star formation at least in 
haloes within this mass range, if not larger, in order to be consistent 
with the observational findings that only a small fraction ( ^ 10 %) 
of baryons is turned into stars at z ^ 3 ([Gawiser et al.|2007[[Moster[ 
et al.|2010[[Guo et al.|2010[[Behroozi, Wechsler & Conroy|2013[ 

Moster, Naab & White|2013[ l. 

The inclusion of thermal SN explosions (NutTFB) has a 
weak impact on the star formation history of our intermediate-size 
galaxy. This is not unexpected, given that thermal energy is suscep¬ 
tible to severe artificial radiative losses. Nevertheless, the galaxy 
stellar mass from the NutTFB run is reduced by 31% (Mgtar = 
1.3 X 10^° M©), compared with that from the cooling run. Fig¬ 
ure (the right panel) shows that the suppression of star formation 
occurs mostly at an early epoch (z ^ 4), suggesting that thermal 
feedback is more efficient in dark matter haloes with a shallower 
potential. 

There is a clear trend that star formation is more suppressed 
in smaller haloes with kinetic feedback (Figures and [^. This is 
essentially because the ejecta momentum from SN explosions with 
a slight boost (p ^ 3pej) is significant enough to disperse small 
star-forming clouds present in dwarf-sized galaxies. However, as 
the halo becomes larger and more gas is accreted, the galaxy forms 
a denser and more massive gas cloud at its centre, leading to effi¬ 
cient star formation. The resulting stellar mass of the Nut galaxy at 
z = 3 in the NutKFB run becomes greater (1.5 x 10^° M©) than 
that in the NutTFB run. 


On the other hand, the use of the appropriate momentum from 
the adiabatic to snowplow phase ensures that star formation is no¬ 
tably suppressed in the NutMFB run. The galaxy stellar mass in 
this run is 2.4 times smaller than that of the NutCO run at z = 3. It 
is worth noting that the mass difference is even more significant at 
high redshifts (z ^ 5), indicating again that the regulation of star 
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Figure 6. The mass of dark matter haloes in which each star particle is 
born. The y-axis indicates the amount of stars formed in haloes of different 
masses. A larger amount of stars form in more massive haloes, consistent 
with the observed main star formation sequence. Star formation should be 
suppressed not only in small haloes but also in massive haloes to be con¬ 
sistent with the empirical results that only a small fraction of baryons turns 
into stars. 


formation is more effective in small-mass haloes. The amount of 
stars formed in haloes with Mhaio ^ M© is reduced by an 

order of magnitude (Figure [^. However, once a giant gas clump 
forms (Mcore ~ 10® M©) as a result of mergers of massive gas 
clumps (Mgas ~ 10^“® M©) at the centre of the halo at z ^ 4.8, 
even momentum from SNe fail to disrupt the cloud (Figure]^. We 
will discuss the cause of the inefficiency in more detail in Section|^ 

In previous cosmological simulations, owing to finite resolu¬ 
tion, assuming that each star particle of mass 10^ — 10® M© ex¬ 
periences a single explosion was unavoidable. It is now feasible 
to model discrete, multiple SNe, at least in zoom-in simulations, 
and to investigate their effect (NutMFBm). There are two impor¬ 
tant differences in MFBm compared to the aforementioned runs. 
First, since SNe start occurring as early as 3 Myr after the forma¬ 
tion of a star particle in the NutMFMm run, the self-regulation pro¬ 
cess starts dispersing dense gas clouds earlier than in the runs with 
a single event (NutTFB, KFB, and MFB). On the other hand, as SN 
explosions are distributed over time, they become less disruptive 
in dense environments. This is precisely the motivation of using a 
“superbubble” feedback in some studies (e.g. |Sharma et al.|2014| l. 
Figuresj^and [^indeed show that the stellar mass grows faster than 
the NutMFB run, indicating that the feedback is not as efficient as 
the single explosion case (NutMFB). 

The last model (NutMFBmp) explores the possibility that SN 
explosions channel preferentially through a lower-density medium 
which may not be fully resolved in our simulations. Note that the 
momentum input from SNe is not boosted artificially in this run, but 
calculated according to Equation with an updated average den¬ 
sity, as a smaller amount of mass from the cell hosting the SN is 
now entrained (/host = 0.1, Equation [^. We find that the subtle 
change in modeling momentum deposition can significantly alter 



z=4.80 (1.273Gyr) z=4.77(1.283Gyr) z=4.71 (1.304Gyr) 


Figure 7. Formation of a massive stellar core due to the overcooling trig¬ 
gered by mergers of gas clumps in the NutMFB run at 2 ; 4.8. The top 

and middle panels display the projected gas density at different redshifts. 
The numbers in the parentheses are the age of the universe at each red- 
shifts. The circles indicate the clumps that lead to a formation of a gi¬ 
ant clump of 1.4 X 10® M© at 2 ; = 4.71 (Mgtar = 8 x 10^ M©, 
Afgas = 6 X 10^ M©) for which cooling is accelerated. The clumps with 
the same colours merge first into a clump before the white and yellow cir¬ 
cles merge at 2 : = 4.77. Also included in the bottom panels are the stellar 
density distributions at 4.71 < 2 ; < 4.80. The escape velocity of the mas¬ 
sive clump is ~ llOkrns”^. Note that winds with v > 110km are 
difficult to generate inside a dense, gas-rich core via momentum input from 
SN. A possible solution to prevent the onset of the overcooling is by having 
early, efficient (SN) feedback. 

star formation histories of the main galaxy from NutMFBmp (Fig- 
ure[^ red lines). As the outflow becomes faster along lower-density 
channels and momentum injection is slightly increased (~309^ if 
the density is reduced by a factor of ten), SN explosions ensure to 
remove gas from dense star-forming regions continuously and sup¬ 
press gas accretion on to the galaxy centre (see Section |^. As a 
result, the galaxy stellar mass is reduced by more than an order of 
magnitude at 2 ; = 3 (Mgtar = 7.4 x 10® M©), compared to that 
of the cooling run. Notice that we are not claiming that the outflow 
along lower-density channels would easily disrupt the massive gas 
clouds with TTigas ^ 10^ M© which already suffer from the over¬ 
cooling problem. Rather, we find that the feedback scheme is able 
to prevent the formation of overcooled structures in the first place 
through the continuous outflows and keep the star formation rate 
low. 

We also note that star formation in the NutMFBmp run be¬ 
comes episodic, as a small number of gas clumps govern the total 
star formation rate (Figure]^. When these clumps are disrupted by 
the effective feedback, star formation drops rapidly. Such bursty 

Note that only a small fraction of cells (10%) would have 30% increase in 
the SN momentum by assuming that the volume-filling density is ten times 
smaller than the mean density of the host cell of SN (see Section]^. The 
total momentum input from SNe in the galaxy is increased only by ~4 %. 
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star formation is also seen in [Hopkins et aL|f2014| Figure 9). At 
lower redshifts (z < 1) when the galaxy becomes more massive 
and its specific star formation rate drops, many star-forming clouds 
are likely to form in the galactic disc (see Figure 1 in [Hopkins et al.| 
[2014[ for example) and star formation histories become smoother. 


It is worth mentioning that the role of the outfiow along lower- 
density channels in regulating star formation may vary under cer¬ 
tain circumstances. For example, using a set of idealised simu¬ 
lations with an initially turbulent medium, [Hennebelle & Iffrigj 
( |2014j ) show that SNe exploding in lower density environments en¬ 
hance star formation (their model D) than the case where SNe al¬ 
ways explode in a dense region (their model Cl, Figure 13). They 
argue that this is likely because most of the SN energy is carried 
away with diffuse gas, compressing the neighboring dense clouds 
(model D). This may sound contradictory to our findings, but we 
note that the different conclusions may be reached due to different 
initial conditions and timescales involved. First, while dense clouds 
develop simultaneously over the entire computational domain in 
[Hennebelle & Iffri^ ( [2014[ ), gas is accreted smoothly through the 
cosmic web over several Gyr in our simulations. Thus, the diffuse 
gas accelerated by SN explosions is less likely to encounter the 
collapsing clouds and trigger star formation in our case. Second, 
although the diffuse hot gas may not destroy nearby star-forming 
clouds, it can reduce the gas cooling and accretion onto the galaxy 
centre by interacting with halo gas, leading to a lower star forma¬ 
tion activity on cosmological time scales. 


In Figurewe examine the stellar mass fraction of the galax¬ 
ies (Mstar//barMhaio) in our simulations. Note that galaxies out¬ 
side the virial radius of the main halo are also included. For com¬ 
parison, we include the empirical determinations of the stellar mass 


fraction based on the abundance matching technique (Moster et al. 

2010|[Guo etal. 

2010|jBehroozi, Wechsler & Conroy[20 13 [[Moster, 

Naab & White 

[2013jl. Independent measurements based on the 


clustering analysis of 3022 star-forming galaxies selected from the 


VIMOS Ultra Deep Surve y ([Durkalec et al.|2015| z ~ 3) and 162 
Ly a emitters at z ~ 3.1 ( [Gawiser et al.|2007[ ) are also included 
as a blue pentagon and a cyan diamond with error bars, respec¬ 
tively. We stress that the comparison between the empirical results 
and our simulations should be taken with caution, given that the 
simulated galaxy sample is quite small. Nevertheless, several inter¬ 
esting features can be gleaned from this plot. First, a large fraction 
of baryons is converted into stars in massive haloes in the cooling 
run, while the fraction is small in haloes of a few times 10® M©, as 
the UV background heating prevents gas from collapsing. We con¬ 
firm that a similar trend is found when the fraction of total bary- 
onic mass fraction is plotted instead of the stellar mass fraction. 
Second, there is a general trend that the stellar mass fraction is 
lower in smaller haloes in the runs with feedback. Third, the use 
of more realistic time delays does reduce stellar mass in the small 
haloes (Mhaio ^ 10^'^ M©). Fourth, the stellar fraction in small 
haloes agrees well between MFBm and MFBmp. This means that 
the fraction of SN host cell mass entrained is irrelevant, as long as 
the momentum injection from SNe is high enough to unbind the 
star-forming cloud. Finally, the MFBmp run produces an amount 
of stars consistent with the empirical sequence, indicating that mo¬ 
mentum input from SNe can provide enough pressure support to 
regulate star formation. Whether this momentum can be efficiently 
transferred to the turbulent medium, as postulated in the MFBmp 
run, remains to be investigated, however. 
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Figure 8. Stellar mass fraction of galaxies in simulations with different 
feedback models at 2 ; = 3. Note that galaxies that are not accreted onto 
the main halo are also included. Different symbols and colours correspond 
to different feedback schemes. Empirical determinations based on the abun¬ 
dance matching technique are shown as different lines, as indicated in the 
legend (Moster, Naab & White[2013[[Behroozi, Wechsler & Conroy|2013). 
The shaded regions denote Icr and 2cr uncertainties for the [Moster, Naab[ 
[& White[ ( [2013) and fBehroozi, Wechsler & Conroy[j20l'3) results, respec¬ 
tively. Independent measurements from the clustering analysis of 2 ; ~ 3 
star-forming galaxies selected from the VIMOS ultra deep survey jPurkalec] 
[et al.|20l3) and Lya emitters at 2 ; ~ 3.1 jGawiser et al.|2007) are shown 
as a blue pentagon and a cyan diamond symbol, respectively. 

It can be seen that star formation in galaxies embedded in lower-mass 
haloes is more suppressed in general. Mechanical SN feedback is more 
efficient than other feedback schemes at regulating star formation. This 
plot demonstrates that momentum input from SNe alone may be able to 
match the empirical sequence (the NutMFBmp run), provided that the ISM 
is highly turbulent. 


3.2 Metal enrichment 


It is well established that the gas phase metallicity in galaxies in 
creases with galaxy lumi nosity or mass ([Garnett[[2002 Pilyugin, 


Vflchez & Contini|2004| [Tremonti et al.|2004[ Kewley & Ellison 


at higher redshifts ([Erb et al.[[2006 [Maiolino et al.[[2008 [Man-[ 

nucci et al.[2009 

[Yuan, Kewley & Richard[2013[ Zahid et al.[2013| 

Wuyts et al.pOl 

4| [Steidel et al.[[2014jl, but a systematic offset is 


found such that more distant galaxies are more metal-poor for a 
given mass. This suggests that metals are not instantly recycled to 
enrich the ISM, but possibly blown away through galactic winds 
( [Dalcanton|2007[ ). 

Stellar metallicity is an important indicator of the overcooling 
problem. [Wise et ^ < [2012a[ ) show that the metallicity of stars is 
very sensitive to the strength of stellar feedback. We also confirm 
that the metal enrichment is directly affected by a galaxy’s ability 
to regulate star formation. In the cooling run, the stellar metallicity 
quickly approaches the solar value (Z = 0.02), as newly synthe¬ 
sized metals are recycled instantaneously within their birth clouds. 
Visual inspection of metallicity distributions within the dark matter 
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Figure 9. Evolution of stellar metallicity in the main galaxy. Different 
colour-codings indicate the runs with different feedback models. Also in¬ 
cluded as black or red dotted lines are the corresponding gas metallicities. 
Stellar metallicity is a good indicator of the overcooling problem, as demon¬ 
strated by [Wise et al.| ( |201^ . We find that the stellar metallicity increases 
rapidly as soon as a dense core forms at the galaxy centre (see also Figure[7j. 
Only the efficient feedback model (MFBmp) forms stellar populations with 
a significantly sub-solar metallicity. 


halo shows that the intergalactic medium (IGM) is barely metal- 
enriched (Figure]^ the top second panel). This demonstrates that 
galaxy interactions are not an efficient process to re-distribute met¬ 
als in this intermediate mass haloes. The galaxies in runs with SN 
feedback start with a lower metallicity when the halo mass is small, 
but again most of them approach solar metallicity as soon as a mas¬ 
sive stellar core which feedback cannot destroy forms in the galaxy 
(see Figure]^. The only exception is the NutMFBmp run in which 
strong outfiows prevent the formation of such a massive core by 
continuously dispersing gas clumps (Figure]^. 

Figure compares the mass-metallicity relation (MZR) of 
the simulated galaxies from different feedback runs at z = 3. We 
compute the gas-phase metallicity by assuming the solar abundance 
((0 /H)q = 8.69, Asplund et al. 2009). Our minimal run (NutCO) 
shows a steep MZR in small galaxies (Mgtar ^10^ M©), as star 
formation is mitigated by the strong UV background radiation. The 
MZR becomes nearly fiat, however, for more massive galaxies, 
again due to efficient recycling of metals, in contradiction to ob¬ 
servations. The MZR has a slope close to the observed one only 
when SN feedback is included. 

The normalization of the MZR in simulations is di¬ 
rectly affected by the strength of SN feedback. While ob¬ 
servations at high redshifts {z ~ 3) show that galaxies 

with 10^ ^ Mstar ^ 10^'^ M© are metal-poor (~ 0.1-0.3Zsun) 
( [Maiolino et al.|200^|Mannucci et al.|2009| l, simulations with weak 
SN feedback (thermal and kinetic) predict metal-rich dwarf galaxy 
populations with ~ 0.3-1.0Zsun. Only the galaxies with star for¬ 
mation regulated by strong SN feedback (e.g. MFBmp) reasonably 
follow the observed sequence. 

It is worth noting that the simulated main galaxies with dif- 


Figure 10. The mass-metallicity relation (MZR) at 2 ; = 3 in the Nut simu¬ 
lations. Different symbols and colours correspond to the runs with different 
feedback models. For comparison, we include the observed MZRs drawn 
from 9 galaxies at z ~ 3.5 jMaiolino et al.|200^ and 8 LBGs at 2 ; ~ 3.1 
JMannucci et al.|2009). The solid line is the fit to the results from|Maiolino| 
leTSIpOSl , while the local MZR derived from the local SDSS galaxies is 
shown as a dashed line jKewley & Ellison|2008) . A fraction of simulated 
galaxies exhibits very enhanced metallicity compared to their mean MZR 
sequence, because SNe blow out star-forming gas temporarily and the gas 
metallicity is governed by the enriched outfiow. Note that without outfiows 
(NutCO) the slope of the MZR is shallower than the observations. The ef¬ 
fectiveness of SN feedback determines the offset from the observed MZR 
at 2 : ~ 3. 


ferent feedback, except the one from the MFBmp run, show sim¬ 
ilar gas-phase metallicities at z = 3. This is essentially because 
the gas metallicity saturates at around the solar value. Although 
the metallicity of the ejecta adopted in this study is 2.5 Zsun, the 
pristine gas is accreted continuously onto the galaxy at high red- 
shifts, diluting the enriched gas. Furthermore, since stellar mass 
increases together with the metallicity (see also |Agertz & Kravtsov| 
|2015| Figure 6), the gas-phase metallicities of the main galax¬ 
ies seem to agree with the observational data within the current 
measurement uncertainties. Given the comparison of stellar mass 
fraction, this means that the gas-phase metallicity of the mas¬ 
sive galaxies (Mstar ^ lO^^M©) may not be a good indicator 
for the overcooling problem, although dwarf galaxies with masses 
10® Mstar ^10^ M© will still be very useful to put a strong 
constraint on galaxy formation models. 

The origin of the MZR is often attributed to the fact that star 
formation efficiency is lower in smaller galaxies. Some authors ar¬ 
gue that outfiows may not be necessary to explain the MZR as long 
as the stellar fraction is a strong function of galaxy mass (|Finlator| 
|& Dave|[2008[ [Mannucci et al.|[200^ [Calura et al.||2009| l. How- 

ever, our simulation (NutCO) suggests that the fraction is always 
high and not a strong function of galaxy (or halo) mass without 
feedback (NutCO, see Figure [^. This is essentially because gas 
experiences a runaway collapse, shifting the density distribution of 
star-forming gas towards a very dense regime. Note that the typical 
density of star-forming clouds without feedback is set by the force 
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Figure 11. Evolution of the half-mass radius (i^eff ,m) of tho main galaxy 
at 3 < 2 ; < 15 (Az ^ 0.5). The galaxy at 2 : = 3 is shown as a bigger 
symbol than others. Different colour-codings and symbols are the same as 
in Figure [To| Galaxies grow with time, as late accretion supplies gas with 
larger specific angular momentum. Strong feedback tends to produce a more 
extended stellar component. Note that the half-mass radii can also increase 
temporarily due to galaxy mergers. The late increase in radius at late times 
in the run with thermal feedback (TFB) is also the outcome of serendipitous 
mergers (see text). 


resolution of numerical simulations, and we find that this is roughly 
nu ~ 10^'^ cm“^ in our runs with 12 pc resolution (see also Fig¬ 
ure 12 of |Hopkins, Quataert & Murray |2012b| for higher resolution 
runs where the typical density is larger nn ^ 10^ cm“^). Thus, 
for a star formation efficiency of Cff = 0.02, one would expect 
that more than 80% of the gas cloud is converted into stars in 100 
Myr. Since the Hubble time is much longer than this for 2 : = 3 
galaxies, metallicity evolves similarly as a closed box model in the 
absence of feedback, except that the “equilibrium” metallicity is 
smaller than the metallicity of the stellar ejecta (Zej = 0.05) due 
to the infall of pristine gas. |Finlator & Dav^ ( |2008| ) already pointed 
out that in the absence of outflows the capability of reproducing the 
MZR depends on the gas consumption time scale (see their Figure 
13). 


3.3 Galaxy size, circular velocity, and morphology 

An important prediction from tidal torque theory is that late accre¬ 
tion carries larger specific angular momentum (j) than the gas in- 
fiow at early times (|Peebles|1969[|Doroshkevich|1970[|White|1984[ 
[Catelan & Theuns|1996| l. Because of t his, the spin parame ter of a 

dark matter halo (A' = jdm/V^RvirVc, Bullock et al.|200l|> can be 
kept nearly constant across different redshifts (|van den Bosch et al. 


2002 [ |Peirani, Mohayaee & de Freitas Pacheco||2004[ Bett et al. 


20071. Numerical simulations indeed confirm that j of dark matter 
near the virial radius is a factor of two to three larger than the net 
j of the halo ( Kimm et al.|2011b||Pichon et al.|2011 [[Stewart et al.| 


|2013[|Danovich et al.|2015|l 

Baryons gain angular momentum by the same large-scale tidal 
torques as dark matter (e.g. |Kimm et al.|2011a) , and thus galactic 
discs are expected to grow by late gas inflows. Figure pT] indeed 
shows that the half-mass radius of the stars (Reff,in) in the main 
galaxy gets larger with time regardless of the feedback model. One 
can also see that the strength of SN feedback generally leads to a 


larger The NutMFB and NutMFBmp runs produce galaxies 

with Reff,in = 0.83 kpc and 1.5 kpc, while the half-mass radii of 
the main galaxy in NutCO, NutKPB, and NutMFBm runs are 0.43, 
0.40, 0.31 kpc, respectively. 

An exception is the thermal feedback run, which produces a 
rather extended stellar component with Reff,ni = 0.93 kpc. We 
find that the relatively large half-mass radius originates from the 
serendipitous mergers at 5 < z < 6. In this run, the merg¬ 
ing satellite galaxies are more massive than those in other runs 
with feedback (Figure [^, and we find that the orbital angular mo¬ 
menta of them are coincidently aligned with the spin of the galac¬ 
tic disc. Consequently, the merger remnant becomes a more ex¬ 
tended clumpy disc (Figure fourth column), even though there 
is still a significant central concentration of mass at its centre. 
However, it is very unlikely that such serendipitous mergers oc¬ 
cur only in certain feedback models. Although we have not made 
dwarves a focus of this study, we confirm that dwarf galaxies 
with 10^ ^ Mstar ^ 10^ M© from the NutTFB run are in fact 
smaller (i^eff ,m ~ 100 pc) than the runs with mechanical feedback 
(Reff,m ~ 500 pc). 

It is also interesting to note that the galaxy size is closely 
linked with the star formation history and metal enrichment in 
stars. The most evident example can be seen from the NutMFB 
run at 4.5 5.5. This corresponds to the epoch during which 

stellar mass and metallicity increase rapidly indicating the on¬ 
set of severe runaway collapse and overcooling (Figure |^. Ac¬ 
cordingly, the galaxy gets significantly smaller during this period 
(8.5 < logMstar/M© < 9.5). 

Another useful galaxy quantity is the circular velocity (Vc(< 
r) = ^GM(< r)/r). Observed bright spiral galaxies are char¬ 
acterised by the smoothly varying baryonic rotation curves (e.g. 
[Kassin, de Jong & Weiner|2006| l. However, Figure shows that 
the maximum circular velocity appears at very small radii (r ~ 
100 pc) in the runs in which the gas-to-stellar conversion efficiency 
of the main galaxy is higher than the empirical estimates (CO, TFB, 


KFB, MFB, and MFBmp), confirming earlier result s ([Navarro & 


Whitel 19931 

|Abadi et al.|2003 Governato et al.|2007||Scannapieco 

et al.|2009| | 

IHummels & Bryan||20I2| |Agertz & Kravtsov||2015|l. 


The peak near the centre sets in earlier and the maximum circular 
velocity is larger in simulations with less effective feedback. The 
velocity curve in the NutMFB run is smoothly rising at z = 7, but 
the peak at the centre begins to emerge at z ^ 5 again due to the 
runaway gas collapse (Figure [^. Only the run with the mechanical 
SN feedback based on the assumption of a porous ISM (NutMF¬ 
Bmp) does not show the peaked circular velocity curve. Note, how¬ 
ever, that the rotation curve of the galaxy turns out to be very dif¬ 
ferent from that of the circular velocity, indicating that there is no 
stellar disc with a well-ordered motion (see Figure|^for the stellar 
density plot in the fourth column). 

In order to see the kinematic properties of the simulated galax¬ 
ies, we quantify the fraction of stars supported by rotation (/rot) in 
Figure We first determine the spin direction of the total stellar 
component within 0.2 virial radius and the rotation velocity (Urot) 
of each star particle along the azimuthal direction. Then the mass 
fraction of stars whose kinetic energy (F^kin,* = mu^/2) is dom¬ 
inated by the rotational component (F^rot,*/F^kin,* > 0.5) is cal¬ 
culated within the inner radius (r < 0.2i?vir), where F^rot,* = 
mUrot / 2. /rot niay be viewed as a proxy for the galaxy morphology 
in the sense that a purely rotating disc would have /rot = 1, while 
a stellar system with only radial motions would have /rot = 0. 
Figure(top panel) shows that the simulated galaxies with a high 
stellar-to-halo mass ratio are marginally supported by rotation. Ap- 
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Figure 12. Circular velocities of the main galaxy from a suite of Nut simu¬ 
lations. The circular velocity is measured as yc(< r) = sJGM{< r)/r, 
where M(< r) is the total mass (baryons+dark matter) inside radius r. Dif¬ 
ferent lines correspond to runs with different feedback models, as indicated 
in the legend. As known in the literature, the runs with weak SN feedback 
show a peaked velocity curve towards the galaxy centre. The run with me¬ 
chanical feedback in a porous ISM (MFBmp) suppresses the central peak 
in the velocity curve. 


Figure 13. Top: fraction of stars supported by rotation (/rot) in the main 
galaxy as a function of the Hubble time. We simply assume that a star is 
rotationally supported if more than half of its total kinetic energy is ac¬ 
counted for by rotational motion. The different colours and symbols corre¬ 
spond to different feedback models, as indicated in the legend. Bottom: /rot 
of the stars younger than 10 Myr. The main galaxy in the runs having more 
stars than the empirical estimates (CO, TFB, KFB, MFB, and MFBm) is 
marginally rotation-supported. Interestingly, we find that the most effective 
feedback (MFBmp) leads to the formation of a spheroidal galaxy. 


proximately 55-70 % of the stars have orbits that are more or less 
aligned {0 < 45°) with the spin axis of the galaxy^iving rise to 
geometrically thick discs (Figure]^ fourth columns 

On the other hand, a large fraction of stars (~ 80%) in the 
main galaxy from MFBmp is supported by random motions, re¬ 
sulting in a spheroidal morphology rather than a disc (Figure]^. 
This seems to be in conflict with the previous claim that strong 
feedback preferentially blows out gas with low angular momen¬ 
tum (e.g. IGovemato et al.|[T007| [Brook et"^|2011||, but c onsis- 
tent with recent high-resolution simulations by |Roskar et al.| ( |2014| ) 
and [Hopkins et al.| ( |2014| l. Using the same grid-based code as ours 
(RAMSES), [Roskar et al.j ( |2014| ) performed idealised simulations 
varying the optical depth to IR photons to investigate the impact of 
radiation pressure on galaxy properties. They find that the model 
with strong stellar feedback, which matches the empirical stellar- 
to-halo mass relation, leads to the formation of a spheroidal galaxy. 
[Hopkins et al.| ( [2014[ ) perform cosmological simulations to study 
the formation of disc galaxies with a sophisticated stellar feedback 
model, including radiation pressure and SN feedback. They show 


^ Note that we do not reorient the galaxy to the edge-on configuration in 
Figure]^ to demonstrate that the spin direction of the simulated galaxies 
may be significantly altered by the feedback model. 


that the morphology of the Milky-Way progenitor at z = 3.4 is 
irregular due to strong outflows (see their Figure 1). Somewhat en¬ 
couragingly, despite the absence of the well-ordered disc at z ~ 3, 
they find that smooth gas accretion at later times leads to the forma¬ 
tion of a realistic, extended disc. In this regard, the spheroidal mor¬ 
phology we obtain from the MFBmp run may not be problematic, 
but we will need to run the simulation down to z = 0 to confirm 
this. 

Although a significant fraction (~ 30%) of stars has orbits 
misaligned with the overall spin of the galaxy, the stars younger 
than 10 Myr form out of a well organised, gaseous disc (/rot ^0.8) 
in the CO, TFB, KFB, and MFBm runs. This can also be seen in the 
projected star formation rate density distributions (Figure fifth 
columns). One may wonder why /rot of the whole stellar popula¬ 
tion does not increase notably for 3 < z ^ 6, given that /rot of 
the young component is high (~ 0.9). This can be attributed to two 
factors. First, the young stars form in the core of the galaxy (~ 300 
pc) where there exists a massive bulge of ~ 2-4x10® Mq. Since 
the motion of the young stars is mainly governed by the gravity of 
the bulge, the interaction with the spherical structure is likely to 
isotropise their motions. More importantly, the kinematically cold 
young stars are susceptible to external disturbances, such as harass¬ 
ment or mergers. In order to check this, we follow the evolution of 
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/rot of the young stars (t < 10 Myr) formed at different redshifts 
{z — 4.5—4.8, Az = 0.05) in the NutTFB galaxy during which the 
spin vector of the (young) stellar component is nearly unchanged 
( ^ 15°), while their /rot drops to ~ 0.5. We find that the evolu¬ 
tionary patterns of /rot in the young population with different ages 
are well synchronised, following the merging histories of star clus¬ 
ters and/or galaxies, irrespective of their age. This means that the 
young stars are continuously transformed into a bulge by the dy¬ 
namical disturbances. 

We find that star-forming gas shows irregular morphologies 
when SN feedback is strong (see Figure MFBmp) or when the 
galaxy is small (z ^ 7 in the MFB runs). This is essentially be¬ 
cause the feedback keeps blowing out/away a substantial amount 
of gas before a kinematically cold gaseous disc forms. Let us sup¬ 
pose that there exists a well-defined, clumpy galactic disc in the 
x-y plane. If star formation is very active, it would drive strong 
outfiows along the z direction destroying the star-forming clouds. 
Even though the gas is initially coplanar, the spin direction of the 
returning gas would no longer be aligned with the galactic disc, 
randomising the direction of the spin direction of the subsequently 
collapsing star-forming clouds. This effect is likely to be more sig¬ 
nificant if bursty star formation takes place in a small number of 
gas clouds, as in the MFBmp run (Figure]^. However, if a gas- 
rich disc becomes massive enough to behave like a buffer region in 
which the returning galactic fountain gas is first mixed with the disc 
gas or other fountain gas, star formation in the well-defined disc is 
likely to continue, increasing the fraction of kinematically settled 
discs { [Kassin et al.|2012| l. 
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Figure 15. The mass-loading efficiency of winds in the main galaxy at red- 
shifts 3 < z < 12. The loading factor is computed by dividing the out¬ 
flow rate measured at the virial radius with the star formation rate averaged 
over 10 Myr. The strong SN feedback model (MFBmp) leads to a high 
/l^out/SFR of 5-10 when the galaxy stellar mass is greater than 10® Mq . 
Such a high mass loading is normally required to reproduce the stellar mass 
functions in large-scale hydrodynamic simulations (e.g. |Oppenheim er et al.| 

|Mo) . 


3.4 Galactic outflows and inflows 

We now turn our attention to the properties of gas outside the main 
galaxy. The inflow or outflow rate is measured by computing the 
flux at radius r as, 

/ITinf(r’) = AtVT Pgas(^j f^)'?^rad(^5 ^)€)( Urad)? 

Mout(r-) = J (iQ47rr^pgas(r;Q)urad(r’;f^)0(urad), 

where Q is the solid angle, Urad is the radial velocity with positive 
meaning outflows, and © is the Heaviside step function. Since we 
are mainly interested in the effects of SN explosions by the central 
galaxy, we exclude the contribution from the satellites by ignoring 
the gas within their tidal radius. 

We find that the cosmic gas inflow is impervious to SN ex¬ 
plosions (Figurep^ panel (a)). Although strong SN feedback (MF¬ 
Bmp) suppresses star formation by an order of magnitude, it can¬ 
not completely destroy the dense gas filaments around the virial 
radius (see Figure]^ first columns), and both the inflow rate and 
the net infalling velocity are only mildly reduced (panel (d)). Other 
runs with weak feedback show a comparable inflow rate as in the 
NutCO run. This suggests that, if a significant fraction of baryons is 
missing from the host dark matter halo of a Milky Way-like galaxy 
( [McGaugh et al.||2010| c.f. |Werk et al.||2014| |, it is unlikely to be 
caused by the suppression of smooth gas accretion, but rather by 
stellar feedback blowing gas out of the halo. 

The more effective a feedback model is at reducing star for¬ 
mation, the larger the outflow rate it produces (panel (b)). While 
the runs with thermal and kinetic feedback (TFB and KFB) blow 
gas out of the virial sphere at a rate of 1-2 M© yr“^, the mechani¬ 
cal feedback runs show more significant outflows of 1-5 M© yr“^ 
(MFB and MFBm) or 5-10 M© yr“^ (MFBmp). The outflow rates 


measured at smaller radii (0.5 i^vir) tend to be larger than those at 
the virial radii (panel (c)), indicating that the outflow is launched 
more or less ballistically. Although we do not include the fig¬ 
ure in this paper, we find that a large fraction (~ 60% ) of 
baryons escapes from the halo in the MFBmp run, resulting in 
(Mstar + Mgas)/Mvir = 0.071 at z = 3 (scc Table [^. Note 
that this is a factor of three lower than that in the cooling run 
((Mstar + Mgas)/Mvir = 0.238f] 


Due to the enhanced outflow and reduced star forma¬ 
tion in the MFBmp run, a high mass-loading factor (77 = 

Mout/Mstar) of 5-10 is obscrvcd for galaxy stellar masses of 
10® Mstar ^ 10^ M© (Figure 751. This is roughly two orders 
of magnitude larger than those of weaker feedback runs (77 ~ 0 . 1 , 
TFB or KFB). Here we use the outflow rate at the virial radius, 
and thus the mass-loading factor could be even higher if it is mea¬ 
sured at smaller radii. We note that such a high wind efficiency in 
dwarf-sized galaxies is required in large-scale cosmological sim¬ 
ulations to match the stellar mass functions ( |Oppenheimer et ah] 


® [Faucher-Giguere, Keres & M^j201 1} also found a similar feature in their 
SPH simulation without feedback, and attributed this to the delayed gas in¬ 
fall due to photoionization by the UV background. However, we find that 
the main halo in the NutCO run shows the high fraction even before the UV 
background is instantaneously turned on at z=8.5. We argue that the high 
baryon fraction is due to the fact that baryons and dark matter particles set¬ 
tle into a dark matter halo in a different way (Kimm et al.|2011a) . Whereas 
gas agglomerates at the centre of a halo once accreted, collision-less dark 
matter particles pass through the central region first and populate the out¬ 
skirts or regions beyond the virial radius of the halo depending on their 
angular momentum. Thus, it may not be surprising that the baryon fraction 
is higher than /^ if one measures the fraction in the most baryon-rich region 
(i.e. halo centre). 
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Figure 14. Properties of gas inflow and outflow around the main galaxy in the simulations. Different feedback models are shown as different colour-codings 
and line styles. The top panels (a-c) display the accretion (dMinf /dt) or outflow rate (dMout / dt) measured at the virial or half of the virial radius. The middle 
(d-f) and bottom panels (g-h) show the flux-weighted radial velocity (Vr) and metallicity (^gas) of the inflow and outflow, respectively. Note that we exclude 
the gas within the tidal radius of satellite galaxies, hence the properties represent the smooth component. Also included in the last panel (i) is the flux-weighted 
hydrogen number density of inflow(dashed lines) or outflow (solid lines) at the virial radius or 0.5 i^vir. respectively (see the legend). For comparison, the 
mean baryonic density of a virialised halo 17Spc{z)Q.\y/Q.rn) is shown as a gray line, where pc{z) is the critical density of the universe. We And that the 
gas accretion rate is more or less similar regardless of the feedback model, indicating that the cold gas filaments are impervious to SN feedback (see Figurej^ 
first columns). The outflow is stronger in the runs in which star formation is more suppressed. However, the flux-weighted velocity and number density of the 
outflowing component at 0.5 i?vir is largely indistinguishable among the different feedback models. The outflow seen in the run without feedback (NutCO) is 
not feedback-driven, but induced by mergers. Finally, we And that the inflowing gas is roughly ten times more metal-poor than outflows. 
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Table 2. Simulated properties of the most massive galaxy embedded in a 1.2 x 10^^ M© halo at 2 ; = 3. From left to right, each column corresponds to the 
model name, the total stellar mass of a galaxy inside 0.2 i?vir of a dark matter halo (Mgtar), star formation rate averaged over the last 10 Myr (M^star), the 
gas-to-stellar mass ratio (Mgas/Mgtar) within the virial radius, the stellar-to-halo mass ratio (Mgtar/Mhaio) within the virial radius, the fraction of baryons 
inside a virial radius (Mbar/Mhaio)^ the mean metallicity of stars (^star), the mean metallicity of gas (Zg^s) within 0.2 i^vir, outflow rate measured at 0.5 
i^vir (Mout), and the half-mass radius of the main galaxy. 


Model 

Ms tar 

[Mq] 

-Mstar 

[Mo/yr] 



Mbtr/A^halo 

^star 

[Z©] 

^gas 

[Zq] 

Mout 

[Mo/yr] 

'Teff ,m 

[kpc] 

NutCO 

1.9 X IQio 

6.2 

0.45 

0.147 

0.238 

1.03 

0.65 

0.02 

0.43 

NutTFB 

1.3 X 10i° 

4.3 

0.69 

0.096 

0.184 

0.95 

0.63 

0.50 

0.93 

NutKFB 

1.5 X IQio 

5.7 

0.58 

0.115 

0.195 

0.98 

0.50 

0.03 

0.40 

NutMFB 

7.8 X 10® 

5.5 

1.31 

0.061 

0.144 

0.62 

0.58 

0.67 

0.81 

NutMFBm 

9.8 X 10® 

2.4 

1.14 

0.062 

0.148 

0.82 

0.30 

0.38 

0.31 

NutMFBmp 

7.4 X 10® 

0.6 

10.2 

0.006 

0.071 

0.06 

0.07 

1.20 

1.5 


|2010| ). [Hopkins, Quataert & Murray | ( |2012a| ) also show using ide¬ 
alised simulations with parsec-scale resolutions that the mass load¬ 
ing is ?7 ~ 10 for their She galaxy where the maximum circular 
velocity of the dark matter halo (86 km s“^) is comparable to our 
main galaxy at 3 ^ 2: ^ 6. The comparison is by no means suffi¬ 
cient to demonstrate the success of the models, however, and future 
works will have to directly compare to observable signatures of the 
outfiows, such as low-ionization absorption lines (jRupke, Veilleuxj 
|& Sanders|2005l|Martin|2005l l. 

We also find that the physical properties of the outflowing gas 
are generally similar in all feedback runs. Figure (panels e,f, 
and i) shows the flux-weighted velocity and hydrogen number den¬ 
sity of the outflow in different runs. The galactic winds (panels e,f) 
are faster (150-300 kms“^) than the inflow (panel d), and their 
density at 0.5 i^vir is comparable to the mean baryonic density 
of a virialised halo 178pc(^)^b/^m), regardless of the feed¬ 
back model used (panel i), where pc{z) is the critical density of 
the universe. This suggests that it is the opening angle of the out¬ 
flow that determines the overall outflow rate. To substantiate this, 
we compare the distribution of the outflowing gas between NutTFB 
and NutMFBmp in Figure The plot shows the Healpix map 
( jGorski et al.|2005| l of hydrogen number density on a sphere of 0.5 
i^vir at z = 3.9 when the outflow velocity of the two runs is very 
similar. For clarity, both the infalling gas and the gas belonging to 
satellites are removed. It is evident from the figure that the outflow 
rate in the MFBmp run is larger, as it subtends a larger solid angle 
than the one in the thermal feedback. 

Note that the typical density shown in panel (i) (Figure is 
sensitive to the radius at which it is measured. For example, it can 
increase up to 0.1-1 cm“^ in the MFBmp run at O.Si^vir, as small 
gas clumps that are blown away from the galaxy can contribute to 
the outflows (Figurethe bottom right panel). These gas clumps 
with riH > 0.01 cm^rnove outwards more slowly (~ 40 km s“^) 
than the mean outflow velocity (~ 90kms“^) at 0.3 i?vir, and 
thus we expect that many of these clumps do not escape from the 
halo (see also jCooper et al.|2008| l. On the other hand, the runs with 
weak feedback (NutTFB or NutKFB) cannot form many outflowing 
clumps, as SN explosions cannot transfer the correct amount of 
momentum to the ISM. 

Figure (panel h) shows that the outflowing gas is metal- 
rich, but generally sub-solar at z > 3 in the main halo. This means 
that SN ejecta entrains a considerable amount of gas, which is at 
least 1.5 times the ejecta mass, even in the weak feedback runs 
(TFB and KFB), given that the assumed metallicity of the SN ejecta 
is 2.5 solar. Thus, a small mass-loading factor of 77 ~ 0.1 seen in 
Figure [^arises simply because not all SN explosions drive large- 


scale winds. Because of the significant entrainment, the metallici- 
ties of the outflow are comparable to the gas-phase metallicity of 
the main galaxy (Figure [^. For example, the main galaxy in the 
MFBmp run with a mass loading factor of 77 ~ 10 reveals a signif¬ 
icantly smaller metallicity of ~ O.IZ© than the ejecta metallicity. 
The cooling run displays somewhat complex behaviour, as the out¬ 
flow in this case is driven by mergers and interactions with satellite 
galaxies. 

As accretion is dominated by the cold filamentary gas in our 
main halo at z > 3 ( jKimm et al.|2011a] l, the metallicity of the in¬ 
flow is found to be an order of magnitude smaller (0.005—0.03 Zq) 
than that of the outflow. The different metallicity between the in¬ 
flow and outflow demonstrates that the outflow is not efficiently 
mixed with the cold flows (see Figure second columns). In¬ 
deed, the metallicities turn out to be different even for the MF¬ 
Bmp run where the outflow rate is comparable to the infall rate. 
Recently, jCrighton, Hennawi & Prochaska| ( |20l'3] ) analyse the ab¬ 
sorption features of a 2; = 2.44 system with the background 
quasar at 2: = 2.66, and find a cool (T < 20,000 K) component 
with 0.007-0.015 solar metallicity. Their neutral column density 
(A^h = ]^q19-50±o.i 6 mctallicity are sur¬ 

prisingly similar with the properties of the inflow in our simula¬ 
tions, supporting their claim that it is a direct detection of the cold 
flow around a LBG. 


4 DISCUSSION 

4.1 Overcooling problem and a viable solution 

In the previous section, we have shown that a large amount of gas is 
converted into stars at the centre of a galaxy if the feedback cannot 
regulate star formation effectively. The overcooling problem is a 
well known issue and reported in many earlier works based on ther¬ 
mal SN feedback ( |Katz|199^ or a kinetic feedback model ( jPuboisj 
jet al.|201^ . In this work, we demonstrated that the use of the mo¬ 
mentum available at the end of the adiabatic phase (NutMFB and 
NutMFBm) cannot eliminate the overcooling problem, even though 
it better suppresses star formation in small galaxies. More impor¬ 
tantly, we find that star formation histories of galaxies predicted 
from hydrodynamic simulations are very sensitive to the details of 
how momentum from SNe is distributed (see also ?). 

In order to understand why some feedback models are more 
effective at blowing gas away than others, we compute the prob¬ 
ability distribution functions (PDFs) of the density at which stars 
primarily form (72H,form, top panel) and the density at which SNe 
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Figure 16. Differences in the opening angle of the outflow between the thermal SN feedback (top) and mechanical feedback run with the assumption of a 
highly porous ISM (bottom). The panels display the distribution of hydrogen number density on a sphere of 0.5 i^vir (=13.3 kpc, left) at 2 ; = 3.9 when the 
flux-weighted outflowing velocity and density are very similar in the two runs. The right panels correspond to the density maps at smaller radii (O.Si^vir = 8.0 
kpc). The gas belonging to satellites galaxies or infalling gas is removed for clarity. The plot demonstrates that the outflow in the MFBmp run is stronger than 
the TFB run, because it subtends a larger solid angle than the former. At smaller radii, the NutMFBmp run drives winds composed of many gas clumps with 
riH ~ 0.1 cm“^. 


explode (nH,SN, middle panel) in FigureThe PDFs of nH,form 
are obtained by computing the actual probability (oc Mgas/tff) of 
star formation in the ISM of the main galaxy from ~50 snapshots 
at 3 < z < 10. The sites of SN explosions are chosen as regions 
where young stars with the age 9 < t < 10 Myr are located for the 
TFB, MFB, and CO runs, or by taking into account the realistic SN 
rates based on the metallicity and age of each star in the case of the 
MFBm and MFBmp runs at 3 < z < 10 ( [Leitherer et al.|199^ . 
Also included in the bottom panel is the mass-weighted PDF of 
the ISM, which we simply calculate for gas within a radius nsM 
[= min (0.2i^vir, 2 kpc)]. The figure shows that in the absence of 
SN feedback, the majority of stars are born in a very dense medium 
(^H,form > 10^ cm“^). Bccausc the stars do not undergo violent 
dynamical events (except for mergers), most of the metals are sim¬ 
ply distributed around their birth cloud (Figure]^. 

The inclusion of thermal feedback lowers the density peak of 
the PDF for star-forming clouds to nH,form ~ 10^ cm“^, and al¬ 
lows 10% of SNe to emerge at densities nH,SN < lcm“^. 
However, most of the SNe explode in the dense ISM (nH,SN ~ 
10^ cm“^). Note that the shell formation radius at this density is 
only ^ 2 pc for metal-poor gas (Z = 0.01 Zq , [Thornton et ah] 
|1998| ). Given that the radius should be resolved by at least three 
computational cells to capture the final momentum from SN ( |Kim| 
|& Ostriker||2015| l, it is not surprising that thermal feedback fails 
to regulate star formation in the simulation with 12 pc resolution. 
Only SNe that explode in an ISM with riH < 2 cm“^ may ap¬ 
proximate the Sedov-Taylor solution with 12 pc resolution, but this 
requires a very stringent refinement that a cell should not contain 
more than 110 Mq of gas. 

The kinetic feedback scheme by [Dubois & Teyssierj 
( [2008] ) is originally designed to alleviate the overcooling 
problem by distributing a boosted momentum of psN = 
y^(l + ? 7 t(;) 2 F^sNMej to the neighbouring cells of a SN, where 
? 7 ^ = min(? 7 ^^,niax, Mhost/Mej) and Mhost is the gas mass of the 
host cell of a SN. However, if the adiabatic phase of SN expansion 
is not resolved, then radial momentum cannot build up and so the 
impact of SNe is essentially limited by the initial momentum in¬ 


jected to the ISM. Furthermore, the assumption used in the model 
that the mass loading only comes from the gas in the host cell of a 
SN reduced the estimation of the loading factor (r]w) in a low den¬ 
sity medium or in regions where there are multiple SNe, compared 
to the MFB run where neighboring gas is also taken into account to 
calculate the swept-up mass. As a result, the sites of star formation 
and supernova explosion in the KFB run turn out to be quite similar 
to those of the TFB run. 

The most notable difference in the mechanical feedback runs 
(MFB, MFBm, and MFBmp) from others is that a significant frac¬ 
tion of SNe explodes in low densities of nn < 1 cm“^ (40%, 
26%, and 78%, respectively). The PDFs of nH,SN are bimodal in 
these runs, which refiect the fact that the low-density gas (nn ~ 
cm“^) becomes more volume-filling in the main galaxy. This 
low-density gas has a wide distribution of temperature ranging from 
10^ to a few times 10^ K, indicating that the multiphase ISM is nat¬ 
urally generated with SN explosions ( jCox & Smith] 1974[|McKee 8^ 
|Ostriker|1977| l. |Kim, Ostriker & Kim ( |2013| ) simulate SN-regulated 
star formation by injecting the final momentum and unlike us, they 
find that the ISM is composed of cold and warm (T ^ 10^ K) gas. 
We attribute this difference to the fact that SNe in their simulations 
are assumed to emerge instantaneously as soon as a dense gas cloud 
of nn = 100 cm“^ forms, whereas the time delay of 3-40 Myr al¬ 
lows SNe to explode not only in the dense environment, but also in 
a low-density medium in our simulations. 

Figure [Tt] demonstrates that a realistic time delay for SNe is 
crucial for preventing gas from collapsing to very high densities 
(tt-h > 10^ cm“^) (the MFBm and MFBmp runs). This is essen¬ 
tially because more frequent explosions keep stirring dense gas. 
The bottom panel indicates that, in the case of MFBm, the mass 
fraction for gas with densities in the range 100 ^ nn ^ 10^ cm“^ 
increases, lowering to nn ^ 300 cm“^ the typical density at which 
a SN explodes. However, the gas mass contained in a 12 pc cell 
with a density of 100 5500 M©) already exceeds the 

mass swept up by the end of the adiabatic phase 1300 M©, for 
0.1 Zq). Thus, the velocity gained from an individual SN would be 
only ~ 2 km if the final momentum of 2.4 x 10^ Mq km 
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Figure 17. Local environments at which star formation occurs (top) or SNe 
explode (middle). The PDFs of the density at which stars form (nH,form) 
are obtained by computing Mgas/% of the ISM in the main galaxy at 
3 < 2 ; < 10. Similarly, we compute the PDFs of the density at which 
SNe explode (nH,SN) using potential SN particles in the main galaxy (see 
text). Note that although there is no energy input, the NutCO run also has 
SNe that release metals at 10 Myr after the birth of a star particle. Also 
included in the bottom panel is the mass-weighted PDF of the gas within 
a small radius [r = min (0.2i^vir, 2 kpc)]. The realistic time delay of 
SNe commencing at 3 Myr suppresses the formation of very dense gas 
(riH ^ 10"^ cm“^, MFBm and MFBmp), allowing SNe to explode at lower 
density (nn ^ 300 cm“^) than the weak feedback runs (TFB or KFB). 
More SNe emerge in lower density environments when SN feedback be¬ 
comes stronger. 


(Equation is equally distributed to 48 neighboring cells of the 
same density. Given that the local escape velocity (~ 8 km s“^) is 
greater than this, it is very unlikely that a single SN would disrupt 
the entire cloud. One way to unbind the gas is for collective winds 
to continually lift up the gas. This may be plausible because the gas 
experiencing a runaway collapse will always produce young stars. 
On the other hand, the ram pressure due to gas accretion onto a 
collapsing gas cloud is likely to counterbalance the expanding mo¬ 
tion. The MFBm run suggests that such collective events do not 
occur frequently enough to unbind the gas. We note, however, that 


this may simply be due to the fact that the structure of the ISM in 
our simulation is intrinsically smooth at our resolution limits (12 
- 48 pc), and that fast outflows along low-density channels are not 
resolved. 

Small-scale turbulence simulations show that the structure of 
the ISM is complex (e.g.|Padoan, Jo nes & N ordlun d|1997[|Ostriker, 
Stone & Gammie|2001[|Molina et al.|2012"nFederrath & Klessen 

2012 | l, and that the volume-fllling gas density is smaller than the 
mean density by exp [—0.5 In (l + + /3))], where /3 

is the ratio of thermal to magnetic pressure, b is the turbulence 
forcing parameter (0.3 < b < 1), and M is the Mach number. 
As aforementioned, it is likely that our simulations underestimate 
the extent to which SNe channel through a lower density medium 
on scales of 12-48 pc ( |Iffrig & Hennebelle|2015||Kim & Ostriker| 
|2015| ), i.e. our resolution scale, hence, too much gas is likely to be 
entrained along with the SN ejecta from the host cell. Reducing the 
mass loading from the host cell, as in the MFBmp run, increases the 
outflow velocity in low-density neighbours, and the input momen¬ 
tum during the snowplough phase increases slightly (oc 

Figure illustrates an example of how much we expect the 
outflow velocity of each cell to be boosted when we assume that 
/w,host = 0.1 (i.e. only 10 % of the gas from the host cell is en¬ 
trained) instead of /w,host = 1. We use the ISM structure of the 
main galaxy at z = 3 from the cooling run to compare the velocity 
for the case where /w,host = 0.1 (ufo.i) to that where /w,host = 1 
(u). We assume that young stars with 3 < t < 40 Myr are SNe. 
The top panel shows that Ufo.i can be a factor of ~ 10 larger 
than V. Note that such a signiflcant increase is rare because it can 
only occur around a very sharp edge of a dense clump where there 
is a rapid spatial variation of the density (pnbor/phost ~ 0 . 01 ). 
We still And, however, that 10% of the gas cells that are di¬ 
rectly affected by SN explosions would have a velocity more than 
double their value in the case in which the ISM structure is as¬ 
sumed to be smooth on the resolution scale. This occurs in gas cells 
where Pnbor/phost ~ 0.1. Of the 10%, the fraction of cells accel¬ 
erated more than the typical escape velocity of the local GMCs (i.e. 
10 km s“^) increases from 25% to 65% for 0.1 ^ nu ^ 10cm“^. 
Thus, there is a higher probability that the gas is blown away 
through the lower density channels. On the other hand, in dense re¬ 
gions with riH ^ 100 cm“^, the velocity change is not signiflcant 
in general, as the local ISM structure is more or less smooth in the 
NutCO run and the mass contained in the neighbour is already sig¬ 
niflcant, compared to that entrained from the host cell of SN. Even 
though only a small fraction of gas gains a substantial velocity, the 
comparison between MFBm and MFBmp indicates that this plays 
an important role in dispersing the star-forming gas (Figure p^. 
We do not And any sign of overcooling in the MFBmp run. This 
leads us to conclude that the propagation ofSN momentum through 
lower-density channels may be crucial in modelling star formation 
in galaxies, providing a viable solution to the overcooling problem. 


4.2 Comparison to previous work 

To single out the impact from SNe, our simulations do not include 
any form of radiation pressure. However, other feedback processes 
come into play before the first SN (t < 3 Myr), and the role of 
SNe in this work may be somewhat exaggerated. Photoionisation 
heating can reduce the density by an order of magnitude around 
young stars by over-pressurising the medium ( [Krumholz & Tan| 
|2007|[SSes et al.|2014| l. Although the determination of the IR op¬ 
tical depth is not straightforward ( [Krumholz & Thompson|[20T^ 
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Figure 18. Expected increase in the outflow velocity of gas when the ISM 
is assumed to be highly porous. The top panel shows the ratio between the 
outflow velocity estimated assuming that all the gas in the host cell of the 
SN is entrained (u) and the velocity in the case where the volume-fllling 
density is ten times lower than the mean density (ufo.i). Approximately 
10% of the cells that are directly affected by the SN explosions would have 
a twice faster outflow along low-density channels than the flducial case. The 
bottom panel represents the velocities of the cells with Ufo.i > 2v. Black 
and orange points denote v and Ufo.i in physical units, respectively. The 
knee present around nn ~ 1 cm“^ corresponds to the transition between 
the adiabatic and snowplough phase. 


[Davis et ar]|2014| |Rosdahl & Teyssier||2015| l, radiation pressure 
can also drive additional turbulence in optically thick regions. In¬ 
cluding the additional feedback mechanisms may further reduce 
star formation than the model with SN feedback alone (e.g. jAgert^ 
jet al.|201^ or the importance of SNe may simply be mitigated, as 
the early processes substitute their role. Bearing this in mind, we 
compare our results with previous studies in this section. 

Using an AMR code, Enzo ( [Bryan et al.||2014] l, [Hummel^ 
[& Bryan] \20\2\ performed 16 cosmological simulations of Milky 
Way-like galaxies with a few hundreds of pc resolution. They in¬ 
cluded SN feedback by heating the host cell of SN. Because the 
thermal energy is quickly radiated away, more than 50% of the 
baryons ends up turning into stars, leading to a too peaked rotation 
curve (see also [Joung, Cen & Bryan|200^ . The resulting star for¬ 
mation history is found to be smooth, in agreement with the results 
from the NutTFB run. Although there are several differences in nu¬ 
merical methods (e.g. hydrodynamic solver) between the two stud¬ 
ies, this suggests that adopting a higher resolution does not solve 
the angular momentum catastrophe ( [Navarro & Steinmetz[[2000[ 
[Governato et al.[[2004[ l so long as the spatial resolution is larger 
than the shell formation radius. 


To compensate for the finite resolution (~50-300 parsecs) of 
cosmological simulations, several studies adopted the cooling sup- 


pression model ([Stinson et al.[2006|[Governato et al.[2007 

[Brooks 

et al.[2011|[Governato et al.[2010|[Avila-Reese et al.|2011 

[Guedes 

et al.|2011| Teyssier et al.|2013|[Shen et al.|2014[). As gas is pre- 


vented from collapsing onto the centre of dwarf to Milky way-sized 


galaxies by the over-pressure from SNe, galaxies in these studies 
showed a smoothly rising rotation curve in general, as in obser¬ 
vations. Because they only presented the results at low redshifts 
{z < 1), it is difficult to make a direct comparison with our work, 
but star formation histories in the Milky Wa y-sized galaxies ap- 
pear rather smooth at high redshifts (z ^ 2) ( Guedes et al.||2011 
c.f. [Shen et al.[2014| for more episodic SF in dwarf galaxies), com¬ 
pared to our MFBmp run or [Hopkins et al.[ \2014\ . [Hummels &[ 
[Bryan[ ( |2012| ) also examined the model in which the radiative cool¬ 
ing is turned off for 50 Myr inalO^^M© halo, and found that the 
central peak in the rotation curve is reduced, but still too many stars 
formed (~ 60%), compared to the local Mstar-Mhaio relation (e.g. 
[Moster et al.|2010[ l. 

Recently, [Hopkins, Quataert & Murray] ( [2011[ [2012a[ l inves¬ 
tigated the effects of momentum transfer from the UV-IR photons 
on the ISM in four different types of isolated galaxies (SMC-like 
dwarf. Milky Way-like spiral, gas-rich star-forming galaxies hosted 
by a ~ 10^^ M© halo, and star-bursting disc in a 10^^ M© halo) 
with an entropy-conserving SPH code ( [Springel et al.|2005] l. They 
carried out a series of simulations, varying the feedback model, 
with parsec scale resolution, and concluded that the radiation pres¬ 
sure is mainly responsible for the generation of galactic winds in 
the actively star-forming, gas-rich galaxies, whereas SN explosions 
play a major role in relatively quiescent galaxies (SMC and MW- 
like galaxies). However, their SN feedback model is based on the 
initial energy of 10^^ erg, which is likely to have underestimated 
its impact. Our Nut galaxy is similar to their gas-rich, star-forming 
galaxy (Sbc), but we find that the SN explosions are capable of 
driving a strong outfiow if the ISM is assumed to be porous (MF¬ 
Bmp). As mentioned, the wind efficiency (Mout / Mstar) from the 
MFBmp run (Figure is comparable or slightly larger than the 
star-forming galaxy (S bc) in[Hopkins, Quataert & Murray [ ( |20 1 2a[ ). 
In a subsequent paper, [Hopkins et al. ( 20I4| | incorporate the final 
momentum from the SN blast wave in a dense medium, in a very 
similar spirit to our mechanical feedback scheme, and concluded 
that both SNe and radiation pressure are essential for reproducing 
the observed relations, such as Mstar-Mhaio relation or Kennicutt- 
Schmidt law ( [Kennicutt[1998[[Bigiel et al.[2008[ l. 

By including SN, stellar winds, and momentum transfer from 
radiation in the Ramses code ( [Teyssier|2002[ l, [Agertz et al.[ ( [2013[ ) 
systematically examined their relative importance in an isolated 
cloud of 10® M© and a disc galaxy embedded in a halo of mass 
10^^ M©. Because of their finite resolution, they used the the lo¬ 
cal size and mass relation of cluster/clumps to estimate the sur¬ 
face density, which is distinct from [Hopkins, Quataert & Murray 


poll I, who directly measured it from their simulations. [Age^ 
[et al.[ ( 2013[ ) pointed out that the inclusion of early feedback is im¬ 
portant for disrupting the dense gas in star-forming regions. Ap¬ 
plying the radiation pressure model to cosmological simulations, 
[Agertz & Kravtsov[ p015[ ) argue that reproducing a realistic disc 
galaxy requires two conditions: i) radiation pressure should be in¬ 
cluded on top of the energy variable that dissipates away slowly 
on a 10 Myr timescale, ii) star formation should be modelled, such 
that a large fraction (e = 10%) of gas is converted into stars on 
a local dynamical timescale. They find that both the Mstar-Mhaio 
relation and the mass-metallicity relation are well accounted for 
if these conditions are met simultaneously. It is interesting to note 
that when a low star formation efficiency (e = 1%) is used, their 
radiation pressure model could not slow down star formation in 
the Milky Way-like galaxy and converted 50% of the baryons into 
stars. The need for a high e appears to suggest that star forma¬ 
tion should be bursty for radiation pressure to control the growth 


© 0000 RAS, MNRAS 000, 000-000 






































































20 T. Kimm et al. 


of galaxies. Given that our MFBmp run achieved such bursty SFH 
with a relatively low star formation efficiency (e = 2 %), the re¬ 
quirement for a high e may be model-dependent. 

Several authors also studied the role of radiation pressure, and 
come to a different conclusion. Using the ART code ( [Kravtsov^ 
[Klypin & Khokhlov||1997| l, |Ceverino et al.| ( |2014| ) argue that ra¬ 
diation pressure from the ionising photons alone provide enough 
energy to unbind star-forming gas clumps, making the role of mul¬ 
tiply scattered IR photons to be less significant. However, the re¬ 
sults are not straightforward to compare with other studies, because 
they model the radiation pressure as a non-thermal term that does 
not dissipate away until the density becomes lower than the star 


formation threshold (nn = 1 cm '^). Agertz et al. (20131 tested a 


similar model by using a “feedback energy variable”, and showed 
that the resulting star formation histories are quite sensitive to the 
dissipation time scale as well as the fraction of the energy that is de¬ 
posited into the non-thermal term. By contrast, pkumer et ^ ( |2013| ) 
claim that a high IR optical depth of ~ 20 is requisite at z > 2 
to form a realistic disc galaxy in SPH simulations. They could rea¬ 
sonably match not only the Mgtar-Mhaio relation but also the size 
and metallicity evolution in 16 dwarf to intermediate-size galaxies 


resimulated from the Aquarius project ( 

Springel et al.||2008j) and 

some simulations presented injOser et al. 

(|2010|). The different con- 


elusions on the role of radiation pressure by UV to IR photons call 
for further tests based on hydrodynamics calculations fully coupled 
with radiative transfer ([Krumholz & Thompson|2012[ |Davis et ^ 

|2014[|Rosdahl & Teyssier|2015| l. 


5 CONCLUSIONS 

We have investigated the impact of SN feedback on the evolu¬ 
tion of the progenitor of a Milky Way-like galaxy, embedded in a 
10^^ M© halo at high redshifts {z > 3). For this purpose, we have 
carried out a suite of zoom-in cosmological hydrodynamic simula¬ 
tion with high spatial (12 pc) and mass (mdm = 5.5 x 10“^ Mo, 
^star = 610 M©) resolution using the adaptive mesh refinement 
code, Ramses, varying the SN feedback model. Our main results 
can be summarised as follows. 

(i) Our simulations confirm that too many stars would form 
when the feedback from SN cannot drive strong outflows. The cool¬ 
ing run without SN feedback (NutCO) turns 85% of baryons into 
stars inside 0.2 i^vir of the main halo (Mhaio ~ 10 ^^ M©), re¬ 
sulting in a much smaller gas-to-stellar mass ratio (Figure than 
observations at high redshifts. The run with thermal feedback (Nut- 
TFB) also shows the overcooling problem, as the cooling radius of 
the SN blast wave is still under-resolved even in our zoom-in simu¬ 
lations. Inclusion of the SN explosion in a kinetic form with a small 
mass loading of rjw < 10 (NutKFB) does not greatly help either, 
as shock-heated gas is again subject to artificial radiative losses and 
the subsequent generation of radial momentum during the adiabatic 
phase is not guaranteed. Nevertheless, the increase in the input mo¬ 
mentum (^1 + r)w) used in the kinetic feedback scheme helps sup¬ 
press star formation in low-mass haloes (Mhaio ^ 10^° M©) by an 
order of magnitude, compared to the cooling run (Figurej^. 

(ii) The mechanical feedback models (NutMFB, MFBm, and 
MFBmp) that transfers the correct amount of radial momentum 
from the adiabatic to snowplough phase is able to substantially 
suppress star formation in small haloes (Mhaio ^ 10^° M©). The 
galaxy stellar mass of the haloes is reduced by roughly two orders 
of magnitude from that found in the cooling run (Figure [^. 


(iii) The prediction of star formation histories in the main galaxy 
depends on the details of how momentum from SN is distributed. 
The model in which SNe emerge relatively late (10 Myr, MFB) 
could not prevent a runaway collapse of a massive cloud at the 
galaxy centre during frequent mergers of clumps (at z 5, Fig¬ 
ure 1^. This led to a rapid increase in galaxy stellar mass from 
2 : ~ 5. The final stellar mass at 2 : = 3 was over-predicted, com¬ 
pared with the typical observational estimates from the clustering 
analysis of LAEs or the abundance matching technique. 

(iv) Inclusion of a realistic time delay for SNe, commencing 
at 3 Myr, efficiently suppresses the build-up of very dense gas 
(nn ^ 10^cm“^) in the main galaxy (FigurejlTl MFBm). But be¬ 
cause gas outflow velocities become smaller, distributing SN explo¬ 
sions in time makes a galaxy more susceptible to the overcooling 
problem than the case where at leastlO SNe explode simultane¬ 
ously (MFB). As a result, star formation is not strongly regulated 
from high redshifts {z ^ 5), and follows that of the weak feedback 
models (TFB or KFB, Figure [^. 

(v) The observed relations, such as the mass-metallicity relation 

and stellar mass to halo mass ratio, are reproduced by the SN feed¬ 
back model with the realistic time delay, provided that the ISM is 
highly porous (MFBmp). The latter assumption is motivated by the 
observation that galaxies at high redshifts are more turbulent than 
their local counterpart ( [Forster Schreiber et al.|2009) . If the ISM is 
highly turbulent, the volume-filling density of the medium would 
be systemically smaller than the mean. This means that taking the 
uniform gas density at the computational grid scale (as in MFBm) is 
likely to underestimate the wind velocities along low-density chan¬ 
nels by entraining too much gas from the host cell of SN. We find 
that reducing the volume-filling density to 10 % (equivalent to the 
supersonically turbulent ISM with the Mach number of ~ 10) suf¬ 
fices to prevent the overcooling issue, yielding a mass loading that 
is ten times larger than the star formation rate (77 10, FigurefTS]), 

as postulated in the momentum-driven winds ( jOppenheimer et ah] 
|2010j ). Although the degree of turbulence is arbitrarily chosen, this 
demonstrates that the momentum input from SNe alone may be 
able to regulate star formation, as observed. However, future works 
based on resolved SN feedback will be required to draw a firm con¬ 
clusion on the relative role of different stellar processes on the evo¬ 
lution of galaxies though. 

(vi) The MFBmp model shows stronger outfiow than the weak 
feedback models (e.g. TFB, Figure [T4](b) and (c)). The difference 
originates from the opening angle of the outflow on the sky (Fig¬ 
ure [T 6 | , not the density or velocity, which are found to be insen¬ 
sitive to the feedback prescriptions (Figure [T4](e), (f), and (i)). On 
the other hand, metallicity of the outfiow is lower in the model with 
stronger feedback and generally follows that of the galactic gas. 
The fact that outflow metallicity (0.07 ^ Z/Zq < 0.7) is much 
lower than the metallicity of the ejecta (2.5Z©) indicates that the 
outflow is efficiently mixed with the ISM before it is blown away. 
The strong outflow in this MFBmp model leads to the formation of 
a spheroid with Mgtar = 7.8 x 10® M© rather than a well-ordered 
stellar disc at 2 ; = 3 (Figure [T^. 

(vii) The cold filamentary inflow around the Nut halo (Mhaio ~ 
10^^ M©) is impervious to SN feedback. The accretion rate at i^vir 
is found to be more or less similar, irrespective of the strength of 
the feedback (Figure [T4| (a)). The metallicity of the inflow (Z = 
0.005 — 0.02 Zq) is an order of magnitude smaller than that of the 
outflow, and not very sensitive to the feedback model. 
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